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I.  Summary 

Previous  semi-annual  technical  reports  under  this  contract  have  dealt 
with  the  generation  of  acoustic-gravity  waves  and  Rayleigh  waves  by 
explosive  sources.  In  particular,  research  was  done  and  reported  in  the 
following  areas: 

a.  Theoretical  acoustic-gravity  waves  generated  by  surface  displacements 
associated  with  models  of  contained  underground  explosions 

b.  Theoretical  acoustic-gravity  waves  and  Rayleigh  waves  generated  by 
underground  explosions 

c.  Preliminary  theoretical  determination  of  M  versus  yield  of  under- 

s 

ground  explosions 

d.  Theoretical  and  observed  distance  corrections  for  Rayleigh  wave 
magnitudes 

e.  Seismic  source  descriptions  of  underground  explosions  and  a  depth 
discriminate 

f.  Radiation  and  attenuation  of  Rayleigh  waves  from  the  southeastern 
Missouri  earthquake  of  October  21,  1965. 

g.  Theoretical  effect  of  yield  and  burst  height  of  atmospheric  explosions 
on  Rayleigh  wave  amplitudes 

h.  Development  of  the  stochastic  inverse  as  applied  to  static  dislocation 
problems 

i.  A  static  dislocation  model  of  1964  Alaska  earthquake 

Section  II  of  the  final  report  is  a  review  d  evaluation  of  source  and 
linear  propagation  theory  of  acoustic-gravity  waves.  The  discussion  is 


limited  to  the  literature  involving  acoustic-gravity  waves  which  travel  in 
the  lower  atmosphere  and  to  solutions  not  involving  finite  difference  and 
finite  element  techniques. 

In  Section  III,  the  static  displacement  data  obtained  for  the  1971 
San  Fernando,  California,  earthquake  are  used  to  try  to  determine  the 
faulting  mechanism  for  that  event.  The  forward  problem  is  solved  using  a 
three-dimensional  dislocation  theory  applicable  to  a  homogeneous  half-space. 
The  dislocation  model  calculated  to  optimally  fit  the  observed  data  indicate 
a  large  degree  of  spatial  variability  of  dislocation  along  the  assumed  fault 
surfaces.  Part  of  the  data  set  used  in  the  inversion  is  that  obtained  by 
converting  changes  in  the  gravity  field  into  changes  in  elevation  associated 
with  the  earthquake. 


SOURCE  AND  LINEAR  PROPAGATION  THEORY 
OF  ACOUSTIC -GRAVITY  WAVES 


This  discussion  is  limited  to  the  literature  involving  acoustic- 
gravity  pressure  waves  which  travel  in  the  lower  atmosphere  and  to 
solutions  not  involving  finite  difference  techniques  such  as  Green 
and  Whitaker  (1968),  Potemra  (1965),  etc. 


One  of  the  earliest  contriuutions  to  the  theoretical  analysis 
of  infrasonic  wave  trains  produced  by  ground  level  explosions  in  an 
atmosphere  was  that  by  Pekeris  (1948).  This  was  closely  followed 
by  Scorer  (1950).  Pekeris  (1948)  formulated  the  problem  of  a  surface 
source  which  injects  vertically  a  volume  of  gas  into  the  atmosphere. 

He  obtained  solutions  for  two  atmospheric  models  in  terms  of  a  Fourier 
and  a  Fourier-Bessel  integrals.  In  one  model  the  atmosphere  was 
isothermal  and  in  the  other  the  atmosphere  consisted  of  a  troposphere 
having  a  constant  lapse  rate,  and  an  isothermal  stratosphere  with  con¬ 
tinuity  of  temperature  at  the  tropopause.  The  roots  of  the  resulting 
dispersion  relations  were  discussed  but  no  attempt  was  made  to  calculate 
the  predicted  wave  train  Scorer  (1950)  using  a  slightly  different, 
but  equivalent,  formulation  calculated  the  pressure  pulse  at  various 
ranges  for  the  two  layer  model  of  Pekeris.  As  pointed  out  by  Hunt, 
Palmer  and  Penny  (1960),  Scorer  sued  a  modified  stationary  phase 
technique  valid  for  very  large  distances  which  was  unlikely  to  give  a 
good  representation  of  experimental  observations. 
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Yamamoto  (1957)  calculated  the  fundamental  mode  of  two-layer  and 
four-layer  atmospheres.  Yamamoto's  calculations  were  motivated  by  his 
observation  that  the  pulses  from  large  explosions  begin  with  long  period 
waves  with  the  later  a. rivals  in  the  dispersed  signal  having  periods  as 
short  as  0.2  min.  Since  Yamamoto  could  not  obtain  a  cutoff  period  of 
less  than  1.5  minutes  for  a  two  layer  model  with  the  upper  layer  colder 
than  the  lower,  he  concluded  that  the  observed  microbarograms  could  not 
be  explained  with  the  two-layer.  For  a  four-layer  model  with  the  top 
layer  being  the  warmest,  he  found  fundamental  mode  solutions  with  periods 
down  to  0.5  minutes  and  a  group  velocity  curve  with  a  minimum.  His  study 
did  not  include  the  pulse  calculations. 

Hunt,  et  al  (1960),  using  the  Pekeris  and  Scorer  model  and  a  two 
isothermal  layer  and  a  three  isothermal  layer  model,  investigated  the 
character  of  the  acoustic-gravity  modes  and  synthesized  their  portion 
of  the  pressure  pulse  produced  by  explosions  on  the  ground.  Their 
discussion  of  the  influence  of  vertical  temperature  structure  on  the 
character  of  these  modes  is  probably  the  most  physically  instructive 
analysis  in  the  literature.  This  is  true  in  spite  of  the  fact  that  in 
adopting  Scorer's  formulation  to  the  isothermal  layer  and  models,  they 
used  the  incorrect  boundary  condition  between  layers  of  continuity  of 
the  perturbed  pressure  evaluated  at  l  he  ambient  position  of  the  interface 
instead  of  at  the  perturbed  position  (Hines,  1960).  To  the  first  order, 
the  two  boundary  conditions  are  equivalent  if  there  is  no  discontinuity 
in  ambient  density  across  the  interface  as  in  Pekeris  and  Scorer.  This 
error  only  effects  Hunt  et  al's  quantitative  results.  Their  qualitative 
conclusions  are  still  valid. 
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For  the  two  isothermal  layer  atmosphere  Hunt  et_  jQ  found  that 
as  long  as  the  stratosphere  half-space  was  colder  than  troposphere 
layers,  there  was  only  one  mode  which  terminated  with  a  high  frequency 
cutoff.  The  velocity  of  this  mode  varied  between  the  two  layer  velo¬ 
cities  and  was  equal  to  the  stratosphere  velocity  at  cutoff,  This  was 
also  predicted  by  Pekeris  and  Scorer.  When  the  temperature  of  both 
layers  is  the  same  and  the  model  represents  an  isothermal  half-space, 
there  is  again  only  one  mode.  This  is  the  Lamb  mode  which  is  trapped  at 
the  rigid  surl  '.e  with  a  vertical  particle  velocity  equal  to  zero  every¬ 
where  and  travels  at  the  acoustic  velocity  of  the  half-space.  Pekeris 
obtained  a  ficticious  second  mode  for  the  half— space  which  corresponds 
to  the  surface  wave  for  a  half— space  with  only  a  free— surface  boundary 
at  its  top.  Once  the  top  layer  was  made  warmer,  they  found  an  infinite 
number  of  modes  with  no  high-frequency  cutoff  for  any  of  the  modes  and 
a  low  frequency  cutoff  for  all  the  modes  except  the  first.  The  analyses 
of  the  roots  of  this  simple  model  demonstrated  that  if  there  is  a  tem¬ 
perature  minimum,  i.e.  sound  channel,  in  the  atmosphere,  then  there  will 
be  no  high  frequency  cutoff.  This  corrected  the  erroneous  "conjecture 
or  Pekeris  (1948)  that  all  stratified  atmospheres  having  a  temperature 
minimum  would  exhibit  a  high-frequency  cutoff.  "  This  model  also  shows 
that  the  higher  modes  will  have  a  low  frequency  cutoff. 

Weston  (1961)  using  a  Green's  function  approach  obtained  the 
solution  for  a  pressure  pulse  produce..  by  explosions  not  only  on  the 
ground,  but  at  various  heights  in  the  atmosphere  of  spherical  earth. 

In  the  formulation  he  shows  how  the  eigenfunction  for  a  plane  stratified 


model  can  be  corrected  for  a  spherical  cap  spreading  to  represent  the 
pulse  propagating  over  a  spherical  earth.  The  pulse  was  synthesized 
for  three  models,  one  the  Scorer  model  and  the  other  two  were  variations 
of  the  stratosphere  minimum  up  to  a  slightly  colder  than  the  surface, 
half-space  beginning  at  50  km. 

Weston  (1962)  calculated  the  pressure  perturbation  as  in  his 
earlier  paper  using  a  more  realistic  model  of  the  atmosphere  with  two 
temperature  ducts  and  a  large  temperature  gradient  in  the  atmosphere. 

From  his  expression  he  determined  there  should  be  a  set  of  modes  with 
velocities  less  than  the  velocity  of  the  more  extreme  temperature 
minimum  in  the  atmosphere  and  with  periods  always  greater  than  the 
Brunt  cut-off  period  associated  with  this  duct.  These  modes  later 
called  gravity  modes  were  missed  by  Hunt  et^  al^  (1960).  Since  the  cal¬ 
culated  pulses  indicated  that  their  amplitude  decrease  drastically  with 
the  calculated  pulses  indicated  that  their  amplitude  decrease  drastically 
with  increasing  burst  height  which  was  in  direct  conflict  with  later 
investigations,  Pierce,  Posey  and  Iliff  (1971)  concluded  that  there  was 
a  coding  error  in  the  computation  of  Weston's  wave  forms.  Other  evidence 
for  this  is  that  the  fundamental  and  higher  acoustic-gravity  modes  of  his 
calculations  had  a  high— frequency  cut-off,  and  according  to  previous 
and  later  investigations,  these  modes  should  not  have  a  cut-off  for  a 
model  with  a  temperature  duct. 

In  order  to  evaluate  the  radial  or  vertical  eigenfunctions  and 
eigenvalues,  Weston  used  numerical  integration  over  the  vertical 
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distribution  of  temperature.  Press  and  Harkrider  (1962),  Pfeffer  (1962), 
and  Pffeffer  and  Zarichny  (1962,  1963)  formulated  the  eigenvalue  problem 
in  terms  of  a  multilayer  model  of  the  vertical  temperature  structure. 
Yamomoto  (1957)  and  Hunt  et_  al^  also  used  an  isothermal  layer  represen¬ 
tation;  but  in  the  multilayer  approach,  the  layer  solutions  are  combined 
in  layer  matrix  such  that  the  multiplication  of  two  contiguous  layer 
matrices  satisfies  the  boundary  conditions  between  the  layers.  This 
approach  is  particularly  suited  for  a  computer  and  does  not  require 
the  investigators  to  limit  themselves  to  only  a  few  layers  as  in  Yamomoto 
and  Hunt  et  al.  In  the  multilayer  investigations  phase  and  group  velocity 
dispersion  curves  and  vertical  pressure  distributions  were  numerically 
evaluated  in  terms  of  the  different  models  of  atmospheric  structure. 

The  results  of  Press  and  Harkrider  (1962)  and  Pfeffer  and  Zarichny  (1962) 
differed  somewhat  because  the  latter  authors  terminated  their  atmosphere 
model  at  a  lower  altitude.  Later  work  with  a  more  complete  model  by 
Pfiffer  and  Zarichny  (1963)  gives  results  which  agree  with  those  of 
Press  and  Harkrider  and  contain  a  wider  variety  of  atmospheric  models. 

The  multilayer  investigations  found  the  same  set  of  modes,  generally 
referred  to  as  acoustic  modes  which  have  no  high  frequency  cutoff  as 
found  in  the  simple  models  of  Hunt  ^t  al^  and  the  gravity  modes  found  by 
Weston  (1962) .  Press  and  Harkrider  found  that  the  sequence  of  maximums 
or  plateaus  in  group  velocity  formed  by  the  set  of  acoustic  modes  and 
the  first  two  gravity  modes  is  sensitive  to  the  properties  of  the 
lower  50  km  of  the  atmosphere.  Ffeffer  and  Zarichny  (1963)  verified 
this  and  indicated  that  the  sequence  of  plateaus  for  the  higher 
acoustic  modes  approach  the  sound  velocity  of  the  lower  wave  guide 
as  frequency  increases. 


In  addition,  the  phase  and  group  velocities  of  the  acoustic  modes  were 
found  to  approach  the  lowest  velocity  in  the  atmosphere  at  the  high 
frequency  limit.  Press  and  Harkrider  (1962).  In  most  models  this  is 
the  upper  guide  at  85  km;  therefore  the  acoustic  group  velocity  minimums, 
which  are  controlled  by  the  phase-velocity  high  frequency  asymptote,  were 
sensitive  in  these  models  to  the  atmosphere  above  .,0  km. 

While  advances  were  being  made  in  the  evaluation  and  understanding 
of  the  atmospheric  wave  guides,  theoretical  investigations  wer^  being 
done  on  modeling  explosive  sources.  Hunt  et_  al^  (1960),  Weston  (1961, 
1962),  Dikii  (1962),  Pierce  (1963),  Harkrider  (1964),  van  Hulsteyn  (1965, 
I  &  II) ,  all  incorporated  a  point  source  in  their  models  for  calculations 
This  Green's  function  source  was  expressed  in  analytic  form  in  the 
frequency  domain  for  an  isothermal  region  by  Dikii  (1962)  and  Pierce 
(1963).  Using  a  different  notation  for  this  source  Row  (1967)  pointed 
out  that  the  waves  emitted  by  this  source  in  an  isothermal  region  were 
dispersive  and  anisotropic,  i.e.  their  velocity  was  a  function  of 
frequency  and  take-off  angle  from  the  vertical.  Tolstoy  (1963)  made 
the  same  observation  using  a  different  approach. 

Harkrider  (1964) ,  using  an  integral  representation  of  the  frequency 
domain  expression  for  the  Dikii  (1962)  and  Pierce  (1963)  point  source, 
synthesized  barograms  for  comparison  with  observed  waves  from  thermo¬ 
nuclear  explosions.  The  agreement  was  satisfactory  but  showed  some 
disturbing  differences: 

1)  The  suppression  of  high  frequencies  with  increasing  yield  was 
much  more  extreme  than  for  that  observed  expecially  for  yields  above 
10  MT. 
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2)  For  theoretical  barograms  in  the  yield,  W,  range  of  1  to  5  MT, 
the  overpressure  varied  as  W'9.  For  larger  yields  the  amplitude  was 
approximately  proportional  to  W6.  The  latter  differed  considerably 
from  the  expected  direct  proportionality  to  yield,  W1'0. 

3)  In  increasing  the  altitude  from  2  to  17  km  altitude,  the  theo¬ 
retical  barograms  showed  an  increase  in  amplitude  by  about  a  factor  of  6. 

Pierce  (1965)  demonstrated  that  the  pressure  scales  f,r  the  synthetic 
barograms  in  Harkrider  (1964),  were  in  error  due  to  the  omission  of  a 
factor  used  in  scaling  source  functions  from  one  altitude  to  another. 

As  Pierce  pointed  out,  the  corrected  effect  of  increasing  source  altitude 
on  the  amplitudes  of  the  early  arriving  signal  is  slight  at  moderate 
heights  and  eventually  at  great  altitudes  the  signal  pressure  decreases 
with  increasing  source  height  as  one  would  expect.  The  net  decrease 
form  2  to  17  km  altitude  is  on  the  order  of  30%. 

In  addition  to  comparison  with  actual  events  Harkrider  (1964),  used 
the  source  excitation  to  verify  the  effects  of  various  portions  of  the 

wave  guide  as  postulated  by  Press  and  Harkrider  (1962),  and  Pfeffer  and 
Zarichny  (1963) . 

The  models  discussed  above  did  not  include  the  effect  of  wind 
structure.  Weston  and  van  Hulsteyn  (1962),  showed  that  the  linearized 
equations  of  hydrodynamics  are  still  separable  if  the  winds  are  hori¬ 
zontal  and  vary  in  direction  and  magnitude  only  with  altitude.  Prid'iore- 
Brown  (1962),  also  derived  the  generalized  equations  for  sound  propagation 
in  an  atmosphere  with  arbitrary  sound-speed  profile  and  horizontal  wind 
profile.  Weston  and  van  Hulsteyn  presented  the  numerical  effect  of  wind 
on  the  phase  velocity  of  one  mode  propagating  downwind.  Pridmore-Brown 
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included  a  steady  state  point  source. 

P'erce  (1965)  extended  the  theories  of  Weston  and  van  Hulsteyn 
and  of  Pridmore-Brown  to  the  propagation  of  infrasonic  waves  from  an 
idealized  point  source  with  an  arbitrary  time  variation.  Pierce  pointed 
out  that  the  approximations  given  by  Pridmore-Brown  imply  that  the  phase 
and  group  velocity  direction  are  the  same  and  point  radially  from  the 
source.  Pierce  using  a  perturbation  method  and  one  mode  showed  that  in 
general  this  is  not  true.  Pierce  (1967  a  and  b)  combined  his  1965 
formulation  with  the  multilayer  approximation  to  explore  in  greater 
detail  the  effects  of  winds  on  modal  propagation.  As  one  might  expect 
for  pure  upwind  and  downwind  directions  where  the  souice  in  line  with 
assumed  direction  of  the  wind  profile,  the  phase  and  group  velocity  are 
both  in  this  direction  and  that  the  primary  difference  between  upwind 
and  downwind  propagation  on  the  early  arriving  portion  of  the  signal 
would  simply  be  that  waves  tra-.eling  downward  travel  faster.  In  addition 
the  maximum  difference  between  ohase  and  group  velocity  direction  occurs 
for  crosswind  propagation. 

Pierce  and  Posey  (1970)  developed  a  barogram  synthesis  program 
using  the  multilayer  approximation  to  a  temperature  and  wind  stratified 
atmosphere.  In  order  to  facilitate  computation  of  the  effect  of  wind  on 
the  pressure  trains  they  neglect  the  perturbation  between  phase  and  group 
velocity  direction  due  to  crosswinds.  They  compared  the  wave  form  for 
two  wind  models  with  two  non-wind  models  in  which  the  sound  speed  was 
equal  to  sound  speed  of  the  wind  models  increased  or  decreased  by 
the  component  of  the  wind  velocity  in  the  direction  of  propagation.  They 
noted  that,  even  though  winds  were  actually  treated  in  a  much  more 
sophisticated  manner  than  simply  using  this  effective  velocity,  the 
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results  imply  that  the  sophistication  had  only  slight  effect  on  the 
predicted  waveforms.  Thus,  as  long  as  the  wind  speed  in  every  layer 
is  much  less  than  the  speed  of  sound  for  that  layer,  the  major  effect 
of  wind  is  simply  to  change  the  effective  sound  speed  profile. 

The  multilayer  approximation  and  the  application  of  normal  modes  to 
atmosphere  have  received  many  criticisms  and  questions  concerning  their 
validity.  Hines  (1965),  pointed  out  that  the  differential  equations 
derived  for  acoustic-gravity  wave  propagation  in  a  continuously  stratified 
atmosphere  have  coefficients  which  depend  on  derivatives  with  respect  to 
height.  Pierce  (1966),  demonstrated  that  the  approximation  is  rigorously 
valid  in  the  limit  of  small  layer  thicknesses.  Pierce  and  Posey  (1970), 
expressed  the  philosophy  of  the  technique  very  well  in  that  the  method  is 
only  a  numerical  integration  technique.  The  atmosphere  is  not  modeled  by 
a  multilayer  model  at  the  outset,  but  only  use  it  as  a  device  to  evaluate 
numerical  quantities  of  interest.  Any  given  multilayer  atmosphere  would 
most  likely  give  erroneous  results  for  sufficiently  short  wavelength 
events.  Guidelines  for  selecting  the  degree  of  layering  have  been  discussed 
by  Pierce  (1967a),  Friedman  and  Crawford  (1968),  and  Vincent  (1969). 

Another  cause  for  concern  was  the  method  terminating  the  atmosphere 
at  altitude.  Even  though  the  dispersion  curves  in  the  frequency  domain 
change  character  as  the  altitude  increases  and  the  boundary  condition  ia 
changed,  the  barograms  for  a  source  and  receiver  in  the  iower  atmosphere 
are  relatively  unaffected  by  the  atmosphere  conditions  below  110  km 
[Harkrider  and  Wells  (1968),  and  Pierce  and  Posey  0970)].  Also  affected 
by  the  upper  boundary  conditions  is  the  long  period  cutofT  above  which 
ong  period  waves  are  no  longer  trapped  in  an  isothermal  half-space. 
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Since  this  represents  vertical  radiation  out  of  the  atmosphere,  an 
atmosphere  terminated  at  altitude  by  a  free  or  rigid  surface  will  no 
longer  have  this  cutoff  and  waves  of  infinite  period  can  be  trapped  in 
the  atmosphere.  Friedman  (1966),  and  Harki’-der  and  Wells  (1968),  per¬ 
formed  numerical  experiments  which  indicated  that  for  the  extreme  case 
of  infinite  period  for  the  near  earth  surface  mode  traveling  in  an 
atmosphere  terminated  at  altitude  by  a  half-rspace,  the  leakage  of  this 
mode  results  in  a  decay  in  ampitude  of  1/3  after  a  distance  equivalent 
to  almost  eight  encirclements  of  the  earth.  Thus  the  method  of  atmosphere 
termination  is  relatively  unimportant  to  near  surface  generated  and 
recorded  waveforms. 

Two  more  questions  concerning  the  modeling  of  the  atmosphere  were  the 
variation  of  gravity  with  height  [van  Hulsteyn  (1965  I  &  II)]  and  the 
neglect  of  velocity  or  temperature  gradients  in  the  isothermal  layering 
approximation.  It  was  felt  by  some  that  the  neglect  of  the  gravity 
variation  as  in  the  programs  of  Pierce  and  Posey  (1970)  might  be  important. 
This  was  unlikely  though  since  the  results  of  Press  and  Harkrider  (1962), 
and  Harkrider  (1964),  were  obtained  using  in  each  layer  a  value  of 
gravitational  acceleration  corresponding  to  that  value  appropriate  for 
the  layer  midpoint  altitude  above  a  spherical  earth,  and  these  results 
were  in  agreement  with  Pierce. 

Strong  positive  gradients  in  the  altitude  range  from  106  to  175  km 
form  anomalous  regions  in  the  atmosphere  where  the  Brient-Vaisala  period 
is  less  than  the  acoustic  cutoff  period.  No  matter  how  fine  the 
layering  is  made  in  an  attempt  to  model  this  reigon,  it  is  theoretically 
impossible  to  obtain  such  a  region  using  isothermal  layers.  Tolstoy  (1967), 
felt  that  the  presence  of  this  region  might  yield  erroneous  predictions 
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based  on  the  multilayer  model.  Greenfield  and  Harkrider  (1971),  studied 
the  effect  of  neglecting  temperature  gradients  by  replacing  these  gradients 
with  layers  having  constant  gradients  and  comparing  the  results  with 
the  multilayer  calculations.  They  concluded  that  the  more  realistic 
linear  temperature  variation  models  increased  the  plateau  velocity  of 
the  fundamental  gravity  mode  by  about  0.3  per  cent  relative  to  previous 
isothermal  models  and  that  the  greatest  effect  is  due  to  modeling  the 
atmosphere  with  a  linear  temperature  variation  in  the  lower  11  km. 

This  W2s  not  surprising  since  in  the  investigated  period  range  this  mode 
is  the  modal  equivalent  of  the  atmospheric  surface  wave  trapped  near  the 
Earth's  surface.  Comparison  of  constant  gravity  models  with  variable 
gravity  models  indicated  that  the  effect  is  less  than  that  caused  by 
using  linear  temperature  variations  instead  of  isothermal  layers. 

Even  with  the  inclusion  of  these  effects,  the  original  anomalies 
between  observation  and  theoretical  barograms  for  large  nuclear 
explosions  still  existed.  A  major  step  in  resolving  these  differences 
was  made  by  Pierce,  Posey  and  Iliff  (1971).  Using  the  energy  source 
formulation  of  Pierce  (1968)  ,  they  calculated  waveforms  and  compared 
them  with  a  number  of  observed  waveforms  previously  exhibited  by 
Harkrider  (1964) .  Their  source  gave  much  better  agreement  in  the 
suppression  of  high  frequencies  with  increasing  yield  for  yields  above 
10  MT.  Their  results  predict  that  the  early  portion  of  waveforms 
received  on  the  ground  are  nearly  proportional  to  yield  and  Increase 
slowly  with  height  of  burst  up  to  a  height  of  the  order  of  40  km  and 
then  rapidly  decrease  with  further  increasing  height. 

Pierce  (1968),  points  out.  that  all  previous  sources  Including  his 
own  (Pierce,  1963),  were  equivalent  to  mass  injection  sources.  The 
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basic  difference  between  the  mass  injection  and  his  energy  injection 
source  in  the  excitation  of  acoustic-gravity  wave  is  that  the  height 
effect  of  the  energy  source  is  controlled  by  the  parcel  pressure  eigen¬ 
function  and  not  the  pressure  eigenfunction  as  for  a  mass  source.  Thus 
according  to  the  single  composite  mode  analyses  of  Pierce  et  al.  (1971), 
the  mass  source  should  have  approximately  the  same  behavior  with 
increasing  altitude  as  the  energy  source  of  their  paper. 

The  Pierce  et  al.  differences  from  Harkrider  (1964),  were  not  due 
to  the  use  of  the  mass  source  but  in  Harkrider 's  requirement  that  the 
pressure  wave  directly  below  the  source  have  a  specified  size  and  time 
history.  This  resulted  in  unrealistic  scaling  of  the  effect  of  source 
altitude  and  yield  because  of  the  great  difference  in  source  altitude 
and  scale  point  altitude  for  large  yields.  For  the  larger  yields  the 
scale  distance  below  the  source  was  much  greater  than  its  altitude. 

Actually  if  Harkrider  (1964)  had  used  the  scaling  technique  of 
van  Hulsteyn  (1965  I  &  IT),  in  which  the  observed  pressure  pulse  Is  used 
to  calculate  an  equivalent  point  source,  the  only  difference  between  his 
formulation  and  Pierce  and  Posey  (1970)  would  have  been  in  the  vertical 
eigenfunction  mentioned  earlier. 

The  concept  of  composite  modes  used  in  Pierce  e_t  al^.  (1971) 
promises  to  give  physical  insight  and  qualitative  prediction  on  the 
behavior  of  complicated  systems.  Press  and  Harkrider  (1962)  noted  that 
the  phase  velocity  versus  period  curves  for  a  realistic  atmospheric 
profile  tend  to  be  arranged  such  that  a  single  curve  of  nearly  constant 
phase  velocity  extending  over  a  wide  range  of  periods  may  be  constructed 
by  the  connection  of  the  horizontal  segments  of  a  number  of  individual 
modal  dispersion  curves.  This  composite  curve  was  identified  by  many  authors 
as  the  real  atmosphere's  counterpart  of  the  edge  mode  predicted  by 
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Lamb  (1910)  for  the  isothermal  atmosphere  model.  Similarly  Pfeffer  and 
Zarichny  (1963)  reported  that  the  step-like  behavior  of  the  pahse  velocity 
curves  thus  appears  to  be  due  to  the  interference  of  two  families  of 
normal  modes  -  one  family  in  which  the  energy  propagates  in  the  lower 
atmosphere  (the  quais-hor izontal  modes)  and  the  other  in  which  it  travels 
in  the  upper  atmosphere  (the  quasi-vertical  modes).  This  has  been  sup¬ 
ported  by  various  numerical  experiments  (Harkrider  and  Wells,  1968)  as 
well  as  by  examination  of  the  height  profiles  of  associated  with  different 
frequencies  and  modes  (Pffeffer  and  Zarichny,  1963). 

The  phenomena  of  composite  modes  is  especially  evident  in  weakly 
coupled  systems  such  as  the  ocean-atmosphere  system.  Press  and  Harkrider 
(1966)  and  Harkrider  and  Press  (1967)  demonstrated  that  the  classical 
incompressible  ocean  surface  wave  with  long  period  asymptote,  /gl[  ,  and 
short  period  asymptote,  ,  where  H  and  A  are  ocean  depth  and 

horizontal  wave  elgnth  respectively,  is  formed  by  such  a  composite  among 
the  higher  order  gravity  models  of  the  atmosphere. 

Some  implications  of  the  identification  of  the  atmosphere  edge 
mode  were  recently  explored  by  Bretherton  (1969)  and  by  Garrett  (1969a, b). 
Bretherton  showed  how  the  mode's  phase  velocity  may  be  found  to  first 
order  in  deviati  ms  of  the  sound  speed  profile  from  that  of  an  Isothermal 
atmosphere  and  to  first  order  in  the  wind  velocity  while  Garrett  extended 
the  method  to  higher  orders  and  thus  succeeded  in  finding  the  dispersion 
of  the  mode.  The  latter  paper  by  Garrett  (1965)  gives  a  detailed 
discussion  of  the  applications  of  the  theory  based  on  the  edge  mode  to 
the  analysis  of  empirical  waveforms,  as  well  as  of  the  limitations  of 
the  theory.  In  spite  of  the  relative  simplicity  of  the  Bretherton- 
Garrett  model,  it  appears  to  be  adequate  for  the  explanation  of  the 
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gross  propagation  speeds  and  dispersive  characteristics  of  empirical 
waveforms,  although  it  does  not  (at  least  to  the  extent  which  the 
perturbation  series  has  yet  been  carried)  explain  more  subtle  features 
such  as  the  inverse  dispersion  suggested  by  the  data  analysis  of 
Donn  and  Ewing  (1962), 

In  Pierce  and  Posey  (1971),  the  Bretherton-Garret t  model  is  further 
developed,  the  objective  being  a  comprehensive  but  simple  model  for  the¬ 
oretical  synthesis  of  transient  waveforms  recorded  at  large  distances 
from  nuclear  explosions.  The  model  developed  there  would  appear  to 
give  a  relatively  simple  explanation  of  (i)  the  magnitude  of  the 
amplitudes  of  the  first  few  peaks,  (ii)  the  periods  between  the  first 
few  peaks,  (iii)  the  fact  that  records  almost  always  begin  with  a 
positive  pressure  rise,  and  (iv)  the  anomalies  in  the  spatial  distri¬ 
bution  of  recorded  peak  overpressures  that  have  been  reported  by  Wexler 
and  Mass  (1962),  (Item  (iii)  above  is  unfortunately  not  widely  appre¬ 
ciated  due  to  the  fact  that  noise  often  obscures  the  early  part  of 
the  record  and  that  some  of  the  published  records  have  inadvertently 
bean  published  with  the  positive  pressure  axis  extending  downwards. 
Professor  Donn  admitted  to  them  that  this  was  the  case  in  the  paper  by 
Donn  and  Shaw,  (1967).)  The  first  three  items  are  also  explained  by 
the  more  elaborate  multimode  theories  (see,  for  example,  Harkrider,  1964) 
and  they  accordingly  demonstrate  that  the  edge  mode  model  gives  results 
similar  to  the  multimode  calculations  of  Pierce,  Posey  and  Iliff  (1971) 
for  the  first  few  cycles.  The  advantage  of  the  edge  mode  model  is  its 
simplicity,  and  this  allows  an  exploration  of  effects  whose  :onsideration 
may  have  been  prohibitively  difficult  within  the  context  of  existing 
multimode  theories.  Such  effects  include  accumulative  far  field  nonlinear 
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effects,  the  effects  of  terrain  variations  and  the  effects  of  horizontal 
atmospheric  variations.  In  their  paper,  methods  of  examining  such  effects 
of  terrain  variations  and  the  effects  of  horizontal  atmospheric  variations. 
In  their  paper,  methods  of  examining  such  effects  via  the  edge  mode  model 
were  outlined,  although  their  detailed  analysis  was  relegated  to  later 
studies. 
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III.  Static  Dislocation  of  the  1971  San  Fernando  Earthquake 

This  section  is  Chapter  IV  of  Ralph  Wilson  Alewine,  Ill's  Thesis 
(1974)  titled  "Application  of  Linear  Inversion  Theory  Toward  the 
Estimation  of  Seismic  Source  Parameters.  References  and  details  cited 
in  this  section  can  be  found  in  the  Thesis.  This  Thesis  has  been  sub¬ 
mitted  and  successfully  defended  as  partial  fulfillment  of  the  requirements 
of  the  Degree  of  Doctor  of  Philosophy  at  the  California  Institute  of 
Technology,  Pasadena,  California. 
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Chapter  4 

A  Static  Dislocation  Model  of  the 
1971  San  Fernando  Earthquake 

^ ■ 1  Introduction. 

On  February  9,  1971,  a  local  magnitude  6.4  earthquake 
was  initiated  below  the  San  Gabriel  Mountains  north  of  the 
San  Fernando  Valley,  California.  The  surface  ruptures 
from  this  event  occurred  some  13-15  km  to  the  south  of  the 
epicenter  along  the  Santa  Susana-Sierra  Madre  fault  zone 
which  forms  the  frontal  system  of  the  southern  boundary 
of  the  San  Gabriel  Mountains.  The  behavior  and  pattern  of 
the  faulting  that  occurred  during  this  earthquake  exactly 
repeated  and  further  developed  the  older,  already  recog¬ 
nized  geologic  features — faulting  and  uplift — of  the 
southern  boundary  of  the  San  Gabriel  Mountains  which  began 
in  mid- Pleistocene  time  (Oakeshott,  1958).  The  stratig¬ 
raphy  and  geologic  structure  of  this  area  have  been  re¬ 
cently  described  by  Wentworth  and  Yerkes  (1971)  and  by 
Proctor  et  suL._  ( 1972 ) . 

In  this  chapter,  we  will  not  dwell  on  the  geologic 
setting  of  the  area  since  it  is  covered  in  some  detail  by 
the  above  references;  rather,  we  will  propose  a  fault 
model,  the  geometry  of  which  we  think  to  be  most  consistent 
with  all  the  geological  and  geophysical  observations  thar. 
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have  been  made  concerning  this  earthquake.  These  observa¬ 
tions  include  the  focal  mechanism  of  the  main  shock,  the 
mechanism  and  spatial  distribution  of  the  aftershock  se¬ 
quence,  and  the  geologic  character  of  the  surface  faulting. 
As  in  the  previous  chapter,  once  we  have  fixed  the  fault 
model  geometry  we  will  then  use  the  static  displacement 
data  available  to  us  along  a  profile  bisecting  the  surface 
trace  of  the  fault  to  determine  the  static  dislocation 
function  for  our  fault  model  which  is  most  consistent  with 
this  data.  Features  of  this  calculated  static  dislocation 
model  will  then  be  examined  for  its  resolvability,  and 
implications  of  the  model  with  its  associated  stresses  and 
strains  will  then  be  discussed.  Let  us  look  first  at  the 
geological  and  geophysical  observations  that  will  constrain 
our  fault  model  geometry. 

4.2  Fault  Model  Constraints. 

Main  Shock  Focal  Mechanism.  Unlike  the  1964  Alaska  earth¬ 
quake,  there  seems  to  be  no  controversy  of  the  choice  of 
fault  planes  derived  from  the  first  motion  studies  of  the  ' 
1971  San  Fernando  earthquake.  The  natural  ambiguity  in 
the  choice  of  fault  planes  (assuming  a  double  couple  mech¬ 
anism)  is  essentially  removed  by  the  geometric  relation  of 
the  determined  hypocenter  and  the  observed  thrust  fault 
surface  ruptures.  Because  the  earthquake  occurred  in  the 
midst  of  a  dense  seismographic  network,  the  fault  plane 
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solution  is  extremely  well  determined  (Whitcomb,  1971; 

Dil linger  and  Espinosa,  1971;  Canitez  and  Toksoz,  1972). 
This  mechanism  is  primarily  thrust  faulting  on  a  northerly- 
dipping  fault.  Whitcomb  et  al^  (1973)  give  the  following 
fault  plane  parameters:  strike,  N67°(+6°)W;  dip,  52°(_+3°) 
to  NE;  rake,  72° (67°-95° ) . 

Surface  Faulting.  With  the  exception  of  the  dip,  these 
parameters  derived  from  the  faulting  in  the  hypocentra] 
region  agree  well  with  those  determined  by  examining  the 
surface  ruptures.  This  surface  faulting  has  been  described 
in  detail  by  a  number  of  authors  and  agencies,  most  notably 
Kamb  et  al^  (197D,  Palmer  and  Henyey  (1971),  and  Proctor 
—  —  ^972).  The  observed  surface  faulting  can  be 
broken  into  three  main  groups  : 

1.  A  fault  line  extending  approximately  westward 
from  the  mouth  of  Big  Tujunga  Canyon  to  Foothill 
Boulevard,  labeled  the  Tujunga  fault  segment. 
Proctor  et  aJU  (1972)  argue  that  the  name  Lakevicw 
thrust  fault  is  the  appropriate  appelation  for 
this  segment.  Nevertheless,  we  will  refer  to  this 
fault  segment  as  the  Tujunga  segment  so  as  to  be 

consistent,  no  matter  if  incorrect,  v-ith  previous 
work . 

2.  A  complex  group  of  short  length  semi-continuous 
faults  located  O.3-O.9  km  to  the  north  of  and 
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sub-parallel  to  the  Tujunga  fault  segment.  This 
group  will  be  referred  to  later  as  the  splay 
fault  group. 

3.  A  fault  zone  extending  approximately  westward 

from  Pacolma  wash  to  Sylmar.  This  fault  zone  is 
slightly  over  1  km  to  the  north  of  the  Tujunga 
segment  and  is  called  the  Sylmar  segment. 

The  surface  scarps  of  the  Tujunga  segment  extend  more 
or  less  continuously  for  11  Vtm  at  an  average  strike  of 
N70°W.  The  observed  dips  on  the  well  defined  scarps  are 
In  the  range  20-25°  to  the  north,  and  the  motion  Is  shal¬ 
low  thrust  and  left  lateral.  In  the  sub-parallel  splay 
group,  the  scarps  tend  to  parallel  the  sediment  bedding 
which  dips  at  about  65°  to  the  north.  These  faults  have 
been  identified  as  bedding  plane  reverse  faults  with  about 
equal  amounts  of  vertical  and  left  lateral  motion  (Kamb 
et  al . ,  1971).  No  direct  measurement  of  the  dip  was 
possible  on  the  Sylmar  segment  which  extends  about  3  km  in 
an  east-west  direction.  Unlike  the  Tujunga  segment,  the 
surface  scarps  are  less  well  defined  here,  and  the  scarp 
displacements  become  more  diffusely  distributed  (up  to 
50  m)  toward  the  western  end  of  this  branch  (U.  s.  Geologi¬ 
cal  Survey  Staff,  1971). 

Aftershocks .  The  aftershocks  of  this  earthquake  have  been 
studied  Intensely  by  a  number  of  researchers.  Hanks  ct  al. 
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(1971),  Wesson  et  al.  (1971),  and  Allen  et  al.  (1971,  1973) 
have  described  the  location  and  depths  of  these  events,  and 
Whitcomb  (1971,  1973)  has  gone  Into  considerable  detail 
Interpreting  the  mechanisms  and  classifying  these  events. 
The  eplcentral  locations  of  the  aftershocks  outline  an  area 
that  Is  roughly  crescent-shaped.  The  locations  are  clus¬ 
tered  on  the  exterior  of  this  outline  Indicating  that, 
except  at  the  edges  perhaps,  the  stress  was  completely 
relieved  In  the  region  where  the  fault  plane  approached 
the  free  surface.  The  distribution  also  suggests  that  If 
we  take  the  edge  of  the  faulting  plane  to  be  defined  by 
the  outer  limits  of  the  aftershock  distribution,  then  this 
edge  Is  well  defined.  The  hypocenter  of  the  main  shock  Is 
located  at  the  lower,  northernmost  edge  of  this  zone  of 
aftershocks.  Hanks  (1972),  Whitcomb  (1971),  Allen  et  al. 
(1973),  and  Whitcomb  (1973)  suggest  that  the  western  edge 
of  the  aftershock  zone  represents  a  steep  flexure  In  the 
fault  plane.  The  aftershocks  along  this  flexure  are  pre¬ 
dominantly  left  lateral  strike  slip  and  tend  to  occur  at 
deeper  depths  than  on  the  main  fault  plane. 

Allen  et  al .  (1973)  give  evidence  for  the  fault  plane 
to  have  an  Intermediate  dip  which  Is  less  than  that  ob¬ 
tained  from  the  focal  mechanism  of  the  Initial  event  but 
greater  than  the  shallow  surficial  dips.  These  authors 
report  that  the  slip  motion  vectors  for  the  aftershocks 
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located  between  the  hypocenter  of  the  main  event  and  the 
surface  ruptures  had  an  average  plunge  of  36°  toward 
N  20°  E.  Additional  supportive  evidence  for  this  inter¬ 
mediate  dip  angle  for  the  fault  plane  between  the  two 
extremities  is  given  by  the  calculated  hypocentral  depths 
of  the  aftershocks.  Whitcomb  et  al .  (1973)  report  that 
these  events  line  up  in  a  zone  which  has  a  dip  to  the  north 
of  between  35°  to  40°. 

4 . 3  Assumed  Fault  Model. 

We  have  assumed  the  following  fault  system  to  conform 
with  all  the  geological  and  geophysical  properties  that 
have  been  determined  about  this  system  which  was  discussed 
previously.  The  fault  system  is  shown  schematically  in 
Figure  4.1.  The  main  thrust  branch  of  this  system  is  com¬ 
posed  of  21  adjacent  fault  elements.  From  the  surface 
down  to  a  depth  of  7  km,  these  elements  define  a  plane 
which  dips  at  a  constant  angle  of  35°.  From  7  km  to  9  km 
depth,  the  main  thrust  branch  has  a  dip  of  40° ,  and  below 
this  depth,  each  of  the  elements  has  a  dip  of  52°.  The 
lower  edge  of  the  fault  system  is  at  a  depth  of  14.5  km. 

The  splay  fault  system  is  modeled  by  a  series  of  imbricate 
reverse  thrust  fault  planes  intersecting  the  main  fault 
plane  in  the  upper  1  km  block.  These  20  fault  elements 
have  a  constant  dip  of  65°.  We  will  see  later  in  this 
chapter  what  effect  the  addition  of  these  splay  faults  has 


Figure  *1.1.  Perspective  view  of  the  assumed  fault  geometry.  The  solid  circles 
indicate  the  mid-point  depth  (km)  of  the  individual  fault  elements;  the  hypo- 
center  of  the  event  is  indicated  by  the  star.  The  block  offset  on  the  right  is 
an  expanded  view  of  the  upper  1.5  km  to  show  the  details  of  the  splay  fault 
group.  The  inversion  profile  extends  from  left  to  right  along  the  surface. 


-116- 


on  the  calculated  dislocation  function.  Also  included  in 
the  fault  system  is  a  model  of  the  San  Gabriel  fault.  This 
fault  is  perhaps  the  major  geologic  feature  of  the  entire 
epicentral  area.  Wesson  and  Wilson  (1972)  have  postulated 
that  this  fault  was  involved  in  the  seismic  activity  that 
accompanied  this  earthquake.  Kamb  et  al .  (1971)  report, 
however,  that  there  was  no  field  evidence  of  displacement 
occurring  on  this  fault  in  the  aftershock  region.  We  do 
note  that  there  is  a  correlation  of  the  area  of  positive 
gravity  changes  and  the  position  of  this  fault.  This 
suggests  that  if  any  displacement  did  occur  on  this  fault 
during  the  earthquake,  it  occurred  at  depth,  and  the  motion 
was  of  the  normal  faulting  type.  This  fault  is  located 
some  7.2  km  north  of  the  main  surface  breaks  on  the  Tujunga 
segment  and  is  modeled  by  a  plane  of  constant,  62°,  dip 
composed  of  6  fault  elements.  In  this  representation,  the 
San  Gabriel  fault  intersects  the  main  thrust  fault  at  a 
depth  of  9  km.  Thus,  the  total  fault  system  is  composed 
of  47  separate  fault  surface  elements.  Individual  elements 
can  be  eliminated  from  the  system  by  merely  setting  the 
weighting  coefficient  for  that  particular  element  to  zero. 

We  have  chosen  to  take  the  length  of  the  fault  ele¬ 
ments  near  the  surface  to  be  15  km.  We  have  interpreted 
the  clustering  of  the  aftershocks  in  the  epicentral  area 
which  form  the  top  of  the  inverted  "U"  to  be  caused  by  the 
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fault  length  becoming  somewhat  less  as  the  fault  goes  to 
greater  depths.  The  fault  lengths  in  our  assumed  model 
taper  from  15  km  at  the  surface  to  9  km  at  the  hypocenter. 
Table  4.1  gives  the  source  parameters  necessary  for  a 
description  of  this  fault  system.  This  table  has  been 
split  into  two  parts.  Part  a  gives  the  parameters  for  the 
fault  elements  which  comprise  the  main  fault.  Part  b  gives 
the  parameters  for  the  splay  fault  group  and  the  San  Gabriel 
fault  approximation. 

In  the  fault  representation  given  here,  we  have  not 
accounted  for  the  1  km  north-south  offset  between  the 
Tujunga  and  Sylmar  fault  segments,  but  rather  have  con¬ 
sidered  the  Tujunga  segment  as  a  continuous  fault  extending 
linearly  over  the  surface  for  15  km.  Because  this  offset 
occurs  close  to  a  lateral  extremity  of  the  fault  plane  and 
we  intend  to  choose  our  data  for  the  inversion  process  from 
near  the  center  of  the  fault  plane,  this  approximation 
should  be  a  good  one.  We  will  see  later  in  this  chapter 
how  this  offset  can  be  modeled  with  the  fault  system 
geometry  given  here. 

4 . 4  Static  Data. 

Although  unfortunate  from  a  humanitarian  standpoint, 
the  earthquake  propitiously  produced  its  maximum  surface 
deformation  in  an  area  of  exceptionally  good  geodetic  con¬ 
trol.  This  deformation  was  mainly  mapped  by  measuring 
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TABLE  4.1a 


Element  # 

Length 

Dip 

Mid  Pt  Depth 

Width 

Ml 

15 

25° 

0.050 

.24 

M2 

15 

25° 

0.150 

.24 

M3 

15 

30° 

0.275 

.30 

M4 

15 

30° 

0.425 

.30 

M5 

15 

35° 

0.600 

.35 

M6 

15 

35° 

0.800 

.35 

M7 

15 

35° 

1.050 

.35 

M8 

15 

35° 

1.350 

.52 

M9 

15 

35° 

1.650 

•  52 

M10 

15 

35° 

1.950 

•  52 

Mil 

15 

35° 

2 . 300 

.70 

Ml  2 

15 

35° 

2.750 

.87 

Ml  3 

15 

35° 

3.250 

.87 

Ml  4 

15 

35° 

3.750 

.87 

Ml  5 

15 

35° 

4.500 

1.74 

Ml  6 

14 

35° 

5.500 

1.74 

Ml  7 

13 

35° 

6.500 

1.74 

Ml  8 

12 

40° 

8.000 

3.11 

Ml  9 

11 

52° 

10.100 

2.79 

M20 

10 

52° 

11.950 

1.90 

M21 

9 

52° 

13.450 

1.90 

Table  4.1a.  Source  Parameters  for  San  Fernando 
Fault.  Parameters  are  in  km. 


Element  # 


51 

52 

53 

54 

55 

56 

57 

58 

59 

510 

511 

512 

513 

514 

515 

516 

517 

518 

519 

520 


SGI 

SG2 

SG3 

SG4 

SG5 

SG6 


Table  4.1b. 
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TABLE  4.1b 


Length 

Dip 

Mid  Pt. 
Depth 

Width 

15 

60° 

0.25 

0.12 

15 

60° 

0.075 

0.12 

15 

60° 

0.04 

0.12 

15 

6  0° 

0.11 

0.12 

15 

6  0° 

0.18 

0.18 

15 

60° 

0.04 

0.12 

15 

60° 

0.13 

0.12 

15 

60° 

0.22 

0.23 

15 

6  0° 

0.31 

0.22 

15 

60° 

0.04 

0.12 

15 

6o° 

0.13 

0.12 

15 

60° 

0.23 

0.17 

15 

60° 

0.31 

0.17 

15 

60° 

0.44 

0.31 

15 

60° 

0.04 

0.12 

15 

60° 

0.13 

0.12 

15 

60° 

0.24 

0.17 

15 

60° 

0.34 

0.17 

15 

60° 

0.46 

0.23 

15 

60° 

0. 60 

0.28 

15 

62° 

0.75 

0.57 

15 

62° 

1.50 

1.13 

15 

62° 

2.50 

1.13 

15 

62° 

3.75 

1.70 

14 

62° 

5.50 

2.27 

12 

62° 

7.75 

2.83 

Source  Parameters  for  San  Fernando 
Fault.  Parameters  are  In  km. 
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elevation  and  horizontal  line  length  changes  by  comparing 
pre-earthquake  and  post-earthquake  measurements.  These 
data  have  been  compiled  principally  by  the  U.  S.  Geological 
Survey  (Burford  et  al .  ,  1971;  Savage  et  al. ,  1973),  the 
City  of  Los  Angeles,  Bureau  of  Engineering  (personal  com¬ 
munication,  1971),  and  the  Survey  Division,  Department  of 
County  Engineer,  County  of  Los  Angelrs  (1971).  From  these 
sources,  over  100  vertical  displacement  data  points  were 
obtained.  These  points  were  then  plotted  and  contoured  for 
equal  elevation  changes.  This  contour  plot  is  shown  in 
Figure  4.2.  With  the  exception  of  a  few  points  located 
just  to  the  south  of  the  surface  ruptures,  all  the  data 
points  obtained  from  the  geodetic  releveling  show  an  up¬ 
lift.  In  addition  to  the  releveling  data,  we  can  estimate 
changes  in  elevation  by  examining  the  changes  in  the 
gravity  field  caused  by  this  change  in  elevation.  For 
this  earthquake,  we  have  measured  gravity  changes  in  un¬ 
surveyed  areas  where  elevation  changes  would  otherwise  be 
unattainable.  Appendix  2  of  this  thesis  describes  the 
gravity  changes  which  accompanied  this  earthquake  and  the 
relationship  between  the  measured  elevation  changes  and 
the  gravity  changes. 

Figure  4.2  also  shows  a  profile  which  approximately 
bisects  the  area  of  deformation.  In  this  study,  wo  will 
limit  our  inversion  data  set  to  vertical  displacement  data 
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Figure  4.2.  Contour  plot  of  observed  vertical 
static  displacements.  Open  circles  indicate  the 
location  of  the  data  values  used  in  the  inversion 
profile.  Profile  is  indicated  as  the  line  per¬ 
pendicular  to  the  strike  of  the  surface  deformation. 


only.  An  attempt  was  made  to  use  horizontal  shortening  in 
the  direction  perpendicular  to  the  strike  of  the  fault, 
but  the  only  data  which  were  least  contaminated  by  end 
effects  were  relative  line  changes  between  stations  both 
of  which  were  located  very  near  the  surface  ruptures. 

These  data  are  shown  in  Figure  4.3  (eg.  TUJ-6P10).  This 
kind  of  data  is  particularly  ill-suited  for  inclusion  in 
the  inversion  because  if  the  surface  ruptures  are  mis- 
located  by  even  10— 20m,  the  relative  line  change  can  change 
drastically.  Since  the  Tujunga  scarp  has  this  much  local 
variation  in  linearity,  we  considered  these  data  too  un¬ 
stable  to  use  in  the  inversion  scheme.  We  therefore 
limited  the  inversion  data  set  to  vertical  displacements 
only. 

As  in  the  case  of  the  Alaska  earthquake,  we  restricted 
our  data  set  to  only  those  points  which  could  be  projected 
onto  a  profile  without  crossing  a  contour  line.  There  were 
16  releveling  points  near  the  profile  which  met  this  cri¬ 
terion.  These  points  are  shown  in  Figure  4.2.  It  is  seen 
on  this  figure  that  along  the  profile  the  contours  for 
displacements  greater  than  100  cm  are  quite  regular.  By 
regular,  we  mean  that  they  are  orthogonal  to  the  profile. 

On  the  50  cm  contour  line  north  of  the  surface  ruptures 
some  asymmetry  is  apparent.  Since  we  wanted  to  include 
information  over  as  much  of  the  profile  as  possible,  we 
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decided  to  Include  an  estimation  of  the  displacement  here. 
The  data  point  between  the  50  cm  contour  line  and  the  0  cm 
contour  line  to  the  west  of  the  profile  shows  35  cm  of  up¬ 
lift,  and  the  data  point  within  the  50  cm  contour  line  to 
the  east  of  the  profile  indicates  65  cm  of  uplift.  The  50 
cm  contour  line  was  drawn  halfway  between  these  two  points, 
and  a  50  cm  data  value  was  included  into  the  inversion 
data  set  at  the  point  where  the  50  cm  contour  line  inter¬ 
sected  the  profile.  In  making  this  approximation,  we  have 
included  a  larger  variance  to  this  point.  We  will  see 
later  that  a  lateral  variation  of  +  0.5  km  of  the  location 
of  this  particular  data  point  along  the  profile  does  not 
degrade  the  fit  of  our  final  model.  There  are  other  re¬ 
leveling  points  along  this  profile  in  the  epicentral  area 
along  Sand  Canyon  Road.  The  vertical  displacement  values 
of  these  points  range  from  25  cm  to  30  cm  over  a  length  of 
2-3  km.  We  have  chosen  to  include  this  cluster  of  points 
as  a  single  value  of  25  cm  and  to  Increase  the  variance 
associated  with  this  datum.  This  was  done  because 
Ellingwood  and  Williamson  (1971)  report  that  these  points 

might  have  as  much  as  8-11  cm  of  uplift  due  to  secular 
change . 

Pour  more  vertical  displacement  data  points  were  added 
to  the  inversion  data  set  by  converting  gravity  changes  to 
elevation  changes.  Three  of  these  points  are  located  to 
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the  north  of  the  points  shown  In  Figure  ^ . 2 ,  but  still 
south  of  the  epicenter.  One  point  is  located  north  of  the 
Sand  Canyon  Road  releveling  points.  The  exact  locations  of 
these  gravity  points  are  shown  in  Figure  A2 .  ^  given  in 
Appendix  2.  Further  discussion  about  the  gravity  points 
can  be  found  in  that  section  of  this  thesis. 

In  this  investigation,  we  will  assume  that  all  the 
vertical  displacement  occurred  at  the  time  of  the  initial 
earthquake.  Displacements  occurring  after  this  main  event 
were  found  by  a  number  of  surveys  to  be  very  minor  when 
compared  with  the  total  vertical  displacement  (Nason,  1971; 
Lahr  et  al,,  1971;  Burford  et  al. ,  1971;  Savage  et  al. , 
1973;  Sylvester  and  Pollard,  1973).  Any  vertical  afterslip 
that  was  measured  for  a  particular  datum  is  included  in  the 
variance  of  the  displacement  for  that  datum. 

Table  it.  2  gives  the  entire  a0ta  set  that  we  will  use 
in  the  inversion.  The  origin  of  the  profile  is  taken  at 
the  surface  breaks  with  positive  distance  taken  to  be  in 
the  direction  of  N  20°  W.  Also  included  in  this  table  is 
the  estimated  error  in  the  data.  This  error  estimate  is 
taken  to  be  that  estimated  by  the  individual  investigators 
who  measured  the  data. 

A.  5  Calculated  Dislocation  Model. 

Starting  Model.  Since  several  other  investigators  (Savage 
et  al . ,  1973;  Canitez  and  Toksoz,  1972)  have  proposed  a 
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TABLE  4.2 


Distance 

(km) 

Az 

(cm) 

6  (Az ) 
(cm) 

Reference 

Year 

Source 

-0.73 

-1.0 

3.0 

I960 

LABE 

-0.23 

-3.0 

3-0 

I960 

LABE 

0.22 

114.0 

4.0 

1963 

USGS 

0.3^ 

142.0 

4.0 

1963 

USGS 

0.50 

185  •  0 

4.0 

1963 

USGS 

0.55 

192.0 

4.0 

1963 

USGS 

1.26 

228.0 

4.0 

1963 

USGS 

1.48 

207 .  o 

4.0 

I960 

LABE 

1.87 

173-0 

4.0 

1963 

USGS 

2. 40 

149.0 

4.0 

1963 

USGS 

2.66 

122.0 

6.0 

1963 

USGS 

2.78 

135.0 

6.0 

1963 

USGS 

3.42 

98.0 

4.0 

1963 

USGS 

4.82 

71.0 

4.0 

1929 

USGS 

6.52 

50.0 

8.0 

* 

USGS 

8.22 

29.0 

13.0 

**  USGS-GRAV 

9.92 

0.0 

13.0 

1970 

GRAV 

11.50 

-15.0 

20.0 

1970 

GRAV 

17.00 

22 . 0 

17.0 

1964  GRAV-USGS 

20.00 

10.0 

17.0 

1964 

GRAV 

*  Average 

-  65  cm 

(1929)  and  38  cm  (i960) 

##  Average 

along  Sand  Canyon  Road  ( 196 3 

-1968) 

Table  4.2. 

Vertical  uplift 

data  along 

Inversion  Pr 

file. 

Distance 

is  measured 

from  origin 

profile 

,  Az  is 

the  uplift. 

6(Az)  is  the 

estimated  uncertainty  of  the  measurement. 


LABE  =  Los  Angeles  Bureau  of  Engineering, 
USGS  =  U.  S.  Geological  Survey,  and  GRAV  = 
Gravity  measurements  (see  Appendix  2). 
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constant  dislocation  on  some  constant-dipping  fault  plane 
to  be  consistent  with  the  static  data  for  this  earthquake, 
one  such  starting  model  would  be  to  assume  a  constant  dis¬ 
location  throughout  the  assumed  main  fault  system.  A 
simple  numerical  experiment  was  devised  to  test  the  applic¬ 
ability  of  this  constant  dislocation  model  to  the  inver¬ 
sion  data  set.  This  experiment  is  described  as  the 
following. 

Prom  Table  k.2  we  see  that  the  maximum  vertical  dis¬ 
placements  occur  at  about  1  km  to  the  north  of  the  origin 
of  our  fault  system.  (The  origin  of  the  fault  system  is 
taken  as  the  point  where  the  Profile  AA '  intersects  the 
line  formed  when  the  main  fault  surface  intersects  the 
free  surface.)  If  we  normalize  the  observed  displacement 
data  to  this  maximum  value,  we  see  that  there  is  a  rapid 
diminution  with  distance  away  from  the  surface  ruptures. 

One  way  of  affecting  this  rapid  distance-decay  of  the  ver¬ 
tical  displacements  is  to  have  the  fault  plane  dip  at  a 
steep  angle.  In  the  case  of  this  earthquake,  the  dip  of  a 
single  fault  plane  representation  is  limited  by  the  spa¬ 
tial  relationship  of  the  surface  ruptures  and  the  epicen- 
tral  location.  If  we  assume  that  the  hypocenter  marked 
the  lower  bound  of  the  fault  plane,  then  the  maximum 
angle  that  a,  single  plane  could  take  would  be  about  ^5°. 
This  takes  into  consideration  the  uncertainty  in  the 
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locations  of  both  the  epicenter  and  hypocenter.  Another 
factor  that  affects  the  vertical  displacement  at  distance 
away  from  the  surface  ruptures  is  the  depth  to  the  bottom 
of  the  fault  plane.  In  this  problem,  this  is  correlated 
with  the  fault  dip.  The  minimum  value  of  depth  that  this 
fault  plane  could  take  would  be  about  9  km.  If  we  calcu¬ 
lated  the  vertical  displacement  from  a  fault  plane  with 
these  parameters,  we  see  that  the  diminution  is  much 
smaller  than  that  exhibited  by  the  data.  These  results 
after  a  similar  normalization  are  shown  in  Figure  4.4. 

For  a  more  reasonable  geometry  of  a  constant  dip  of  35° 
and  a  depth  of  faulting  of  10  km,  the  fit  to  the  data  is 
even  worse.  Intermediate  to  these  two  curves  in  this 
figure  is  that  for  a  fault  plane  dipping  at  45°  and  extend¬ 
ing  to  a  depth  of  12  km. 

One  might  argue  that  local  geologic  heterogenities 
might  cause  such  a  departure  from  the  idealized  elastic 
case  for  which  we  have  assumed  in  these  calculations.  This 
same  experiment  can  be  carried  out  with  the  finite  element 
modeling  technique.  When  the  geologic  structures  of  this 
area  are  used  in  this  technique,  practically  the  same  re¬ 
sults  are  obtained  (Jungels,  personal  communication,  1972). 

If  we  take  a  planar  dislocation  surface  of  constant 
dip,  but  allow  the  displacement  to  vary  along  this  surface, 
we  can  fit  the  observed  vertical  displacement  much  better 
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if  the  displacement  is  allowed  to  vary  piece-wise  with 

2 

depth  in  a  manner  e  .  No  inversion  was  used  in  these 
calculations,  but  rather  different  slip  models,  each  hav¬ 
ing  a  different  slip  decay-rate  away  from  the  surface,  were 
tried.  The  constant,  y,  necessary  to  fit  the  data  in  this 
case  can  vary  from  0.07  km"2  to  0.12  km"2  depending  on  the 
chosen  dip  of  the  fault  plane.  Mikumo  (1973)  has  arrived 
at  a  similar  result. 

This  simple  calculation  has  been  rather  unsophisti¬ 
cated,  but  we  think  that  this  illustration  clearly  shows 
that  the  observed  vertical  displacement  data  require  a 
fault  dislocation  function  which  varies  along  the  fault. 
This  variation  in  the  fault  slip  should  be  in  the  direc¬ 
tion  of  decreasing  with  increasing  depth  in  order  to 
properly  fit  the  data.  Since  we  have  normalized  the  amp¬ 
litudes  to  that  obtained  about  1  km  distance  from  the  point 
at  which  the  fault  surface  intersects  the  free  surface, 
this  conclusion  might  not  hold  for  the  upper  1  km  (approx¬ 
imately)  of  the  fault  surface. 

We  can  now  proceed  with  the  actual  inversion  of  the 
displacement  data  to  get  the  "best  fit"  dislocation  func¬ 
tion.  We  will  use  as  a  starting  model  one  in  which  the 
displacement  varies  in  an  exponential  manner  as  that 
described  above.  We  can  note  here  that  we  would  Intui¬ 
tively  expect  that  the  displacement  data  shown  in  T'igure  4.4 


would  be  affected  by  only  the  upper  5-7  km  of  the  fault 
system  that  we  have  chosen.  This  fault  system  was  made 
to  consist  of  many  small  fault  elements  in  order  to  more 
adequately  estimate  this  large  slip  variation  which  we 
will  expect. 

Inversion  Results.  The  media  response  matrix, A,  was  cal¬ 
culated  from  the  analytic  expressions  given  by  Mans  inha  and 
Smylie  (1971)  for  the  structural  model  given  in  Table  L\  .1 . 
Using  this  formalism,  the  static  displacement  response  on 
the  surface  at  the  positions  of  the  observed  data  due  to  a 
unit  (1  cm)  dislocation  on  each  of  the  individual  fault 
elements  was  calculated.  Since  there  were  20  data  points, 
or  20  positions  on  the  surface  at  which  the  displacement 
response  was  calculated,  the  media  response  matrix  was 
composed  of  20  rows.  The  number  of  columns  of  the  matrix 
indicates  the  number  of  fault  elements  which  were  used  to 
represent  the  fault  system.  This  number  was  varied,  as 
we  shall  see  below. 

Initially,  we  assumed  that  the  fault  system  was  com¬ 
posed  of  only  the  21  elements  listed  in  Table  H . la  (Mi¬ 
l'll).  We  will  refer  to  this  fault  system  as  the  main 
thrust  sheet.  The  media  response  matrix  for  tills  fault 
system  was  composed  and  the  inversion  scheme  initiated 
with  the  exponental  decay  starting  model  described  in  the 
previous  section.  The  inversion  converged  quickly  to  a 
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"best  fit"  solution.  This  solution  is  presented  in 
Figure  4.5. 

The  lower  part  of  Figure  4.5  shows  this  "best  fit" 
dislocation  on  each  of  the  21  fault  elements  which  com¬ 
prise  this  present  fault  system.  We  see  that  the  disloca¬ 
tion  at  the  surface  is  calculated  to  be  266  cm.  This  value 
immediately  decreases  to  160  cm  at  a  depth  of  about  0.3  km. 
The  maximum  dislocation  (515  cm)  on  the  fault  surface  is 
at  a  depth  of  0.8-0. 9  km.  At  depths  greater  than  1  km,  the 
dislocation  on  the  fault  surface  decays  almost  linearly. 

A  minimum  is  reached  at  about  8 . 5  km  depth,  and  the  dislo¬ 
cation  function  increases  slightly  in  the  hypocentral 
region.  The  maximum  displacement  in  the  hypocentral  region 
in  this  case  is  140  cm.  The  resolvability  of  the  features 
of  this  particular  model  will  be  discussed  in  a  later- 
section  . 

The  displacement  at  the  surface  due  to  this  particular 
dislocation  function  was  calculated  along  a  profile  at  0.2 
km  increments.  This  displacement  profile  is  shown  in  the 
upper  part  of  Figure  4.5.  The  observed  displacement  data 
as  projected  onto  this  profile  are  indicated  along  with 
the  estimate  of  their  errors.  In  this  figure,  we  see  that 
for  the  observed  data  located  to  the  north  of  the  surface 
breaks  (given  by  positive  distances  along  the  abelssa  In 
the  figure),  the  fit  as  given  by  this  "best  fit”  model  Is 
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Figure  tJ.5.  Lower — Best  fit  dislocation  function 
for  the  main  thrust  sheet .  Upper — Comparison  of 
the  calculated  surface  displacements  from  this 
model  to  the  observed  profile  data. 
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quite  good.  To  the  south  of  the  surface  breaks,  the  cal¬ 
culated  model  predicts  too  large  a  negative  surface  dis¬ 
placement  than  Is  observed.  A  noteworthy  feature  of  this 
particular  fit  Is  the  Importance,  as  far  as  spatial  cover¬ 
age  along  the  profile  Is  concerned,  of  having  included  the 
gravity— converted  displacement  points  In  the  Inversion  data 
set.  These  converted  data  points  force  the  calculated  dis¬ 
placement  field  to  go  through  a  slight  minimum  in  the  pro¬ 
file  distance  10-13  km  away  from  the  surface  breaks.  The 
displacements  at  distances  greater  than  at  the  point  where 
this  minimum  occurs  are  dependent  almost  totally  on  the 
dislocation  function  In  the  hypocentral  region.  The  RMS 
error  residual  (observed-calculated)  for  all  the  data  is 
about  8  cm. 

In  the  next  figure  (Figure  4.6),  we  see  what  effect 
there  is  on  this  "best  fit"  solution  when  we  exclude  from 
the  inversion  data  set  the  gravity-converted  displacement 
data.  In  this  figure,  we  see  that  the  dislocation  function 
along  the  upper  6  km  of  the  fault  system  Is  essentially 
unchanged.  Below  this  depth,  the  dislocation  minimum  along 
the  fault  sheet  is  reduced  from  what  was  previously  calcu¬ 
lated  and  shifted  to  a  slightly  smaller  depth.  Larger 
dislocations  (maximum-275  cm)  are  now  given  in  the  hypo- 
central  region.  However,  we  note  that  only  one  data  point, 
located  at  17  km  along  the  profile,  is  controlling  the 
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amount  of  dislocation  in  this  region.  The  resolvability 
of  the  dislocations  in  the  hypocentral  area  will  thus  be 
highly  dependent  on  the  variance  of  this  one  data  point. 

In  the  section  on  averaging  in  this  chapter,  we  will  see 
how  the  inclusion  of  the  gravity-converted  displacement 
data,  inaccurate  as  they  are,  affects  the  resolution  in  the 
hypocentral  area. 

Nex:,  we  want  to  look  at  the  effect  of  including  a 
representation  of  the  San  Gabriel  fault  into  our  fault 
system  for  which  we  will  use  the  inversion  scheme  to  obtain 
a  dislocation  model.  The  structural  parameters  for  this 
addition,  which  is  composed  of  6  elements,  are  given  in 
Table  4.1b  (SG1-SG6).  The  media  response  from  these  ele¬ 
ments  was  computed  and  included  in  the  response  matrix, 
figure  4.7  shows  the  results  of  the  inversion  of  the  data 
set  (including  the  gravity-converted  points)  for  this 
particular  fault  system.  Comparing  the  dislocation  func¬ 
tion  found  in  this  figure  and  that  obtained  by  the  single 
thrust  sheet  (Figure  4.5),  we  see  that  practically  the 
only  difference  between  the  two  models  is  the  amount  of 
fault  dislocation  in  the  hypocentral  area.  In  this  present 
model,  the  slip  in  the  hypocentral  area  is  about  C0% 
greater  than  that  shown  in  Figure  4.5.  This  larger  dis¬ 
placement  in  this  region  is  necessary  to  accommodate  the 
slip  calculated  for  the  San  Gabriel  fault.  The  slip 
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Figure  k.J.  Lower--Best  fit  dislocation  function 
for  the  main  thrust  sheet  and  San  Gabriel  fault 
representation.  Upper — Comparison  of  the  calcu¬ 
lated  surface  displacements  from  this  model  to  the 
observed  profile  data. 
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Figure  ^.8.  Best  fit  solution  for  normal 
the  San  Gabriel  fault  representation. 
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calculated  for  the  elements  representing  the  San  Gabriel 
fault  is  not  shown  in  Figure  4.7  but  is  indicated  in 
Figure  4.8.  The  sense  of  motion  of  the  slip  calculated  on 
the  San  Gabriel  section  is  opposite  that  on  the  main  thrust 
sheet,  i.e.  normal  faulting.  We  see  here  that  the  normal 
slip  decays  rapidly  almost  linearly  toward  the  surface. 

This  model  predicts  no  motion  on  this  fault  for  the  upper 
2  km,  agreeing  with  the  lack  of  observable  expressions  of 
surface  motion  on  this  fault  after  the  earthquake  (Kamb 

1971).  A  later  section  of  this  chapter  goes  into 
detail  about  whether  this  slip  is  required  by  the  data. 

Some  mild  supportive  evidence  for  the  existence  of 
normal  faulting  on  the  San  Gabriel  fault  representation  is 
given  in  the  following  discussion.  Whitcomb  et  al .  (1973) 
give  a  strong  correlation  to  the  noted  areas  of  compres- 
sional  release  (normal  faulting)  as  determined  by  focal 
mechanism  studies  and  the  area  of  mapped  gravity  increases. 
If  the  areas  of  the  San  Gabriel  fault  representation  that 
were  computed  to  have  undergone  normal  faulting  were  pro¬ 
jected  to  the  surface,  then  this  projection  would  approxi¬ 
mately  coincide  with  the  area  of  gravity  Increase  and  the 
correlated  compressional  release. 

Calculation  of  the  averaging  operator  (equation  (2.29)) 
for  this  structural  model  reveals  a  strong  coupling  between 
thrust  slip  in  the  hypocentral  area  and  normal  slip  along 


the  lower  part  of  the  San  Gabriel  fault  representation. 
This  coupling  is  in  the  sense  that  in  order  to  keep  the 
same  displacement  at  the  surface  while  increasing  the 
thrust  displacement  in  the  hypocentral  region,  it  would  be 
necessary  to  increase  the  normal  slip  along  the  lower  part 
of  this  auxiliary  fault.  This  coupling  is  illustrated  by 
the  results  shown  in  Figure  4.9.  In  this  figure,  we  have 
significantly  increased  the  dislocation  in  the  hypocentral 
area.  Hanks  (197*0,  on  the  basis  of  an  analysis  of  the 
Pacoima  Dam  strong  motion  records,  postulated  that  this 
earthquake  initiated  with  large  (4-9  m)  ruptures  in  the 
hypocentral  region.  In  this  example,  we  will  see  If  the 
surface  static  displacement  data  support  this  hypothesis. 

In  this  calculation,  we  have  initially  set  the  fault  dis¬ 
location  for  the  three  bottom  fault  elements  of  the  main 
thrust  sheet  at  the  levels  indicated  in  Figure  4.9.  In 
the  inversion  calculations,  the  model  weights  for  these 
elements  were  set  such  that  these  three  dislocation  values 
remained  fixed.  The  result  from  the  inversion  is  the 
smooth  model  given  in  the  figure.  The  fit  to  the  observed 
data  for  this  particular  model  Is  extremely  good,  with  the 
RMS  error  residual  found  to  be  less  than  5  cm.  We  see  that 
the  minimum  in  the  fault  dislocation  function  is  shifted 
to  a  much  shallower  depth  than  in  previous  examples.  The 
dislocation  function  calculated  for  the  upper'  part  of  the 
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Figure  4.9.  Lower — Best  fit  dislocation  function 
for  the  main  thrust  sheet  and  San  Gabriel  fault 
representation  when  large  displacements  are  forced 
in  the  hypocentral  region.  Upper— Comparison  of 
the  calculated  surface  displacements  from  this 
model  to  the  observed  profile  data. 
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fault  sheet  Is  the  same  as  before.  However,  the  maximum 
normal  slip  on  the  San  Gabriel  fault  has  been  Increased  to 
almost  300  cm  —  approximately  double  that  shown  In  Figure 
4.8.  The  dislocation  function  on  the  auxiliary  fault  re¬ 
tained  the  relative  shape  as  that  shown  In  Figure  4.8. 

We  next  include  the  splay  group  of  fault  elements  into 
our  fault  system.  The  structural  characteristics  of  these 
elements  are  listed  in  Table  lb  (S1-S20).  The  media 
response  to  these  elements  was  calculated,  and  the  response 
matrix  was  expanded  to  include  these  values.  For  compari¬ 
son,  we  chose  as  our  starting  model  for  the  inversion  the 
final  model  given  in  Figure  4.9.  The  results  of  this  In¬ 
version  are  shown  in  Figure  4.10. 

A  comparison  of  the  slip  model  in  this  figure  and  that 
shown  in  the  previous  figure  reveals  that  the  two  models 
are  practically  identical  except  for  the  upper  1. 5-2.0  km. 
This  is  to  be  expected  since  the  fault  geometry  was  only 
changed  for  the  upper  1  km,  and  we  would  expect  this  small 
part  of  the  fault  system  to  be  relatively  uncoupled  from 
the  rest  of  the  fault  system.  For  this  model,  the  fit  to 
the  data  to  the  north  of  the  surface  ruptures  is  equally 
satisfactory  as  with  the  unsplayed  fault  geometry,  but  the 
fit  to  the  two  data  points  to  the  south  of  the  surface 
ruptures  is  degraded  somewhat.  For  the  dislocation  on  the 
main  thrust  sheet  near  the  surface,  instead  of  varying 
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Figure  4.11.  Best  fit  solution  for  slip  on  the 
near-surface  splay  faults. 
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rapidly  as  in  Figure  4.9,  the  slip  is  almost  constant  at 
250  cm,  and  the  maximum  slip  on  the  thrust  sheet  is  re¬ 
duced  by  about  20$  to  400  cm.  The  slips  on  the  small  splay 
faults  vary  smoothly  away  from  the  value  on  the  main  thrust 
sheet  where  they  intersect  that  sheet.  The  slip  on  the 
splay  faults  is  better  illustrated  in  Figure  4.11.  Here 
we  see  that  not  only  does  the  calculated  slip  decay 
smoothly  toward  the  surface,  the  displacements  at  the  sur¬ 
face  decrease  away  from  the  surface  ruptures. 

Since  we  have  limited  our  inversion  data  set  to  points 
which  were  located  near  the  profile  which  approximately 
bisects  the  observed  surface  deformation,  we  would  now  like 
to  see  how  the  calculated  displacement  from  this  model 
compares  to  what  is  observed  in  areas  away  from  this  pro¬ 
file.  For  observation  points  away  from  the  center  of 
deformation,  one  must  consider  the  contribution  to  the 
vertical  displacement  caused  by  any  strike  slip  component 
of  fault  dislocation.  In  this  calculation,  we  will  arbi¬ 
trarily  assume  that  the  dislocation  slip  vector  on  each  of 
the  fault  elements  has  a  constant  rake  angle  of  60°.  This 
is  about  the  maximum  component  of  left  lateral  motion  that 
vvould  be  allowable  based  on  the  focal  mechanism  data  of 
Whitcomb  et  al .  (1973)-  At  this  rake  angle,  the  left 
lateral  strike  slip  dislocation  is  0.6  of  tne  pure  dip 
slip  component.  For  the  majority  of  the  observations,  the 


calculated  vertical  displacement  from  the  strike  slip  com¬ 
ponent  of  dislocation  is  a  secondary  effect. 

One  obvious  feature  of  the  observed  surface  vertical 
displacements  is  the  offset  between  the  Tujunga  and  Sylmar 
fault  segments.  This  offset  is  easily  modeled  by  slightly 
altering  the  geometry  of  the  upper  part  of  the  fault 
system.  An  advantage  of  including  the  small  splay  faults 
near  the  surface,  in  addition  to  approximating  the  surface 
ruptures  observed  north  of  the  main  Tujunga  scarp,  is  that 
it  allows  one  to  easily  model  the  observed  offset  In  the 
surface  faulting.  This  is  done  by  letting  the  upper  six 
elements  of  the  main  thrust  sheet  and  all  the  splay  fault 
elements  (S1-S14),  except  the  ones  fartherest  from  the 
main  surface  breaks,  terminate  at  the  Tuj unga-Sylmar  off¬ 
set.  This  means  that  these  particular  fault  elements  have 
lengths  of  only  xO  xm  with  the  coordinate  system  origin  of 
these  elements  shifted  2.5  km  to  the  east.  The  total  con¬ 
tribution  to  the  surface  deformation  from  all  the  fault 
elements  was  calculated  taking  into  consideration  this 
lateral  offset  where  appropriate.  The  results  are  shown 
in  Figure  4.12.  The  fault  dislocation  function  used  for 
these  calculations  is  the  same  as  that  given  in  Figure  4.5 
with  the  addition  of  the  splay  fault  system. 

A  comparison  of  the  contours  presented  in  this  figure 
and  the  areal  displacement  given  in  Figures  4.2  and  A2.4 
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Figure  4.  12 .  Computed  vertical  displacement  field 
from  best  fit  model  given  in  Figure  ',.5  with  the 
addition  of  the  near-surface  splay  faults.  The 
fault  dislocation  includes  a  component  of  left- 
lateral  strike  slip  taken  to  be  0.6  that  of  the 
dip  slip  motion.  Contours  are  in  cm  of  uplift. 

The  solid  circles  show  the  location  of  the  data  used 
in  the  inversion.  Superposed  on  the  contours  is 
a  map  of  the  roads  in  the  area. 
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shows  che  degree  of  fit  the  calculated  model  has  to  the 
observations.  If  instead  we  had  used  the  dislocation 
function  like  that  indicated  in  Figure  4.7,  the  effect  on 
the  area  deformation  would  be  to  increase  the  magnitude  of 
the  downward  displacement  in  the  northeast  area  and  to  in¬ 
crease  the  northward  extent  of  the  10  cm  contour  line.  For 
contour  levels  25  cm  and  above,  there  is  practically  no 
discernible  difference  between  the  two  plots. 

One  feature  of  Figure  A2.4  which  does  not  appear  in 
the  theoretical  contour  plots  is  the  delineation  of  a  zone 
of  downwarping  to  the  west  of  the  profile  line.  One  poss¬ 
ible  explanation  for  the  occurrence  of  this  area  la  that 
it  is  associated  with  the  seismic  activity  along  the  seis- 
mically  mapped  down  step  of  the  main  fault  plane  (Whitcomb 
et  al- »  1973).  A  significant  amount  of  left  lateral  motion 
on  a  steeply  dipping  plane  could  produce  the  localized 
downward  displacements  that  are  observed  in  this  area.  The 
addition  of  this  fault  geometry  was  not  consilored  in  this 
study,  however. 

4 • 6  Resolvabllity  of  Features  in  the  Slip  Model . 

The  variance  operator,  V,  for  this  case  is  a  47  x  47 
matrix.  This  matrix  was  calculated  in  the  same  manner  in 
the  previous  chapter.  For  simplicity,  the  splay  fault 
elements  were  suppressed  in  these  resolvability  calcula¬ 
tions.  This  suppression  should  cause  no  change  in  the 
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conclusions  reached  in  this  section.  Again,  to  test  the 
resolvability  of  a  particular  feature  of  a  model,  a  model 
perturbation  affecting  this  feature  will  be  considered. 
Tit's  perturbation  will  then  be  tested  for  detectability  at 
the  95%  confidence  level. 

The  initial  model  that  we  consider  is  that  shown  in 
Figure  4.5.  The  first  question  that  we  ask  about  this 
slip  model  is  the  following:  "Considering  that  the  fault 
geometry  of  this  model  is  appropriate,  do  the  data  demand 
that  there  exists  a  sharp  minimum  in  the  slip  in  the  upper 
1  km  of  the  fault  sheet?"  To  answer  this  question,  we 
assume  a  perturbation,  which  when  added  to  the  calculated 
slip  model,  will  just  cancel  this  minimum.  This  perturba¬ 
tion  is  then  tested  and  uniformly  reduced  until  it  is  at 
the  95%  confidence  limit  detectability  threshold.  The 
resulting  perturbation  is  shown  in  Figure  4.13a.  We  see 
that  the  maximum  perturbation  allowable  by  the  data  Ls  not 
large  enough  to  erase  this  minimum. 

Next,  we  try  to  determine  if  the  data  demand  the 
existence  of  the  broad  minimum  in  slip  between  the  upper 
and  lower  parts  of  the  thrust  sheet.  Figure  4.13b  shows 
the  amount  of  slip  that  could  be  added  in  this  region.  We 
see  here  that  the  Silip  gradient  could  be  smaller  than  that 
given  in  the  original  slip  model,  but  nevertheless,  the 


Figure  4,13.  Examples  of  the  resolvability  of 
features  of  the  best  fit  model  solution  given 
in  Figure  4.5. 
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dislocation  does  become  smaller  with  increasing  depth  in 
this  region. 

Intuitively,  we  would  expect  that  the  resolvability 
would  be  less  in  the  hypocentral  area  than  in  other  parts  of 
the  fault  system  because  the  fault  elements  in  this  region 
are  located  farther  from  the  points  on  the  surface  where  we 
have  data.  Also  the  density  of  data  points  on  the  surface 
near  the  epicenter  is  smaller  than  near  the  surface  ruptures 
and  the  data  in  this  area  are  generally  not  as  accurate. 

In  testing  the  resolvability  of  the  calculated  dislocation 
in  the  hypocentral  area,  we  will  consider  what  Is  the  maxi¬ 
mum  perturbation  that  we  can  subtract  from  the  calculated 
values  and  still  be  resolvable  by  the  data.  From  equation 
(2.33)  we  see  that  the  sign  of  the  perturbation  to  be  tested 
by  our  resolvability  criterion  Is  not  important,  only  the 
magnitude  of  the  perturbation.  Thus,  we  can  either  add  or 
subtract  a  perturbation  tested  at  a  particular  confidence 
level  to  the  calculated  model  dislocation  and  still  not  be 
detected  by  the  data  at  that  confidence  level.  In  Figure 
1*.13c,  we  see  that  we  can  add  or  subtract  about  100  cm  of 
displacement  in  the  hypocentral  region  for  this  perturbation 
shape.  From  this  type  of  analysis,  we  indeed  find  that  the 
slip  in  the  hypocenter  Is  poorly  determined,  since  the  un¬ 
certainty  in  slip  is  almost  the  equivalent  of  the  calculated 
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values.  If  required  to  give  bounds  for  the  slip  in  this 
area  as  predicted  by  the  static  data  used  here,  we  would  say 
that  the  maximum  slip  Is  approximately  130  cm  ±  100  cm. 

We  now  consider  the  variance  operator  In  which  the 
San  Gabriel  fault  representation  Is  Included  In  the  fault 
geometry,  but  the  splay  system  Is  still  excluded.  The  slip 
model  for  which  perturbations  are  considered  Is  that  given 
In  Figures  4.7  and  4.8.  Figure  4.14  shows  these  perturba¬ 
tions.  In  the  upper  part  of  this  figure,  we  see  that  con¬ 
sidering  displacement  on  the  San  Gabriel  fault  alone,  that 
Is,  not  taking  Into  account  any  coupling  between  fault 
elements,  there  would  be  about  a  100  cm  uncertainty  In  the 
estimate  of  the  slip.  The  lower  part  of  this  figure  shows 
the  maximum  undetectable  perturbation  possible  when  slip 
on  this  auxiliary  fault  and  slip  In  the  hypocentral  region 
are  considered  simultaneously.  We  see  here  that  the  un¬ 
certainty  In  slip  In  the  hypocentral  area  Is  almost  the 
same  as  that  In  Figure  4.13c,  while  the  uncertainty  on  the 
San  Gabriel  fault  has  been  reduced  to  about  50  cm. 

4 . 7  Averaging  Operators. 

As  discussed  In  the  previous  chapter,  examination  of 
the  kernels  of  the  matrix  operator,  R,  provides  Information 
about  the  filter  through  which  we  are  able  to  view  our 
best  fit"  estimate  of  the  fault  slip.  We  say  that  the 


Figure  4.14.  Examples  of  the 
features  of  the  best  fit  model 
Figures  4.7  and  4.8. 
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ability  to  resolve  the  details  of  the  actual  dislocation 
function  depends  on  the  size  of  the  kernels  and  the  aver¬ 
aging  widths  of  the  kernels.  The  ratio  of  the  size  of  the 
kernels  to  the  averaging  widths  provides  a  convenient 
medium  for  expressing  the  resolution  at  a  given  point  in 
the  fault  system.  First  of  all  in  this  section,  we  would 
like  to  examine  these  ratios  to  determine  the  effect  of 
the  addition  of  the  gravity-converted  displacement  points 
to  the  data  set.  The  averaging  operator  was  calculated 
both  with  and  without  these  data  points  in  the  inversion 
data  set  and  the  averaging  ratios  calculated  in  both  in¬ 
stances  . 

Figure  4.15  shows  the  difference  in  the  averaging 
ratios  the  addition  of  this  data  causes.  The  fault  geometr 
considered  in  this  case  is  that  of  a  single  thrust  sheet. 

In  the  figure,  we  see  that  the  upper  part  of  the  fault  is’ 
unaffected  by  the  addition  of  this  data.  The  area  of  the 
fault  which  has  the  maximum  effect  is  in  the  range  from 
7-11  km  depth.  In  this  range,  the  resolvability  ratio  is 
doubled.  The  reason  for  the  increase  is  not  so  much  an 
increase  in  the  value  of  the  averaging  kernels  for  these 
elements,  but  more  so  the  fact  that  the  averaging  widths  of 
these  elements  are  reduced.  Thus,  the  addition  of  these 
data  does  not  help  much  in  the  determination  of  the  actual 
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Figure  4.15.  Averaging  ratios  for  data  sets  with 
and  without  the  gravity-converted  displacement  data 
Curve  has  been  made  continuous  by  connecting  values 
calculated  at  the  fault  element  mid-point  depths. 


slip  in  this  area,  but  the  slip  that  is  determined  is  much 
more  localized  than  could  be  determined  without  these  data. 

Figure  4.16  shows  examples  of  the  averaging  kernels 
for  the  San  Fernando  earthquake  model  consisting  of  a  sim¬ 
ple  thrust  sheet.  The  arrows  in  the  figure  indicate  the 
elements  around  which  the  averaging  is  centered.  The  block 
offset  to  the  right  is  an  expanded  view  of  the  upper  1.5  km 
of  the  cross  section  to  elucidate  the  fault  geometry  in  that 
area.  Examining  this  figure,  we  see  that  the  averaging  is 
very  localized  in  the  upper  part  of  the  fault,  but  the 
averaging  length  increases  with  increasing  depth  along  the 
fault.  In  the  bottom  figure,  the  averaging  length  (half- 
width)  is  about  4.5  km.  In  comparing  the  averaging  kernels 
presented  here  with  those  given  for  the  Alaska  earthquake, 
we  notice  that  the  strong  negative  side  lobes  are  not  as 
in  evidence  In  this  present  e*  lple.  This  is  probably 
caused  by  the  better  spatial  attribution  of  data  in  the 
San  Fernando  case. 

Figure  4.17  Indicates  the  averaging  kernels  for  se¬ 
lected  fault  elements  when  the  San  Gabriel  fault  geometry 
is  introduced.  The  top  figure  is  centered  on  a  fault  ele¬ 
ment  that  is  located  on  the  main  thrust  sheet  directly 
under  the  San  Gabriel  fault  representation.  We  see  from 
this  picture  that  there  is  no  coupling  in  this  area  between 


Figure  4.17.  Spatial  averaging  of  slip  in  model 
with  San  Gabriel  fault  representation. 
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the  main  fault  and  the  auxiliary  fault.  This  means  that 
there  can  be  no  trade-off  in  displacement  between  the  two 
faults  to  affect  the  resolvability  test  of  the  slip  in  this 
area  as  indicated  in  Figure  4.13b.  Figure  4.17b  shows  the 
spatial  averaging  of  slip  for  the  upper  part  of  the  San 
Gabriel  fault.  We  see  here  that  this  too  is  uncoupled 
from  the  rest  of  the  fault  system.  We  note  here  that  the 
best  fit  solution  for  the  fault  geometry  which  included 
this  fault  indicated  little  or  no  slip  on  this  portion  of 
the  fault.  Figure  4.17c  indicates  the  extent  ^  the  spa¬ 
tial  averaging  for  the  lower  part  of  the  San  Gabriel  fault. 

We  see  here  that  there  is  very  strong  coupling  to  the  fault 
elements  in  the  hypocentral  area.  We  note  that  the  aver¬ 
aging  between  the  two  fault  segments  is  in  the  same  direc¬ 
tion  indicating  that  any  normal  slip  on  this  portion  of  the 
San  Gabriel  fault  will  be  seen  by  the  data  as  norma]  slip 
averaged  over  both  the  San  Gabriel  fault  and  the  main  fault 
in  the  hypocentral  area.  However,  we  can  rule  out  the 
possibility  of  the  occurrence  of  normal  faulting  on  the 
main  fault  from  consideration  of  the  focal  mechanism  of  the 
initial  motion.  These  observations  are  taken  to  mean  that 
this  portion  of  the  San  Gabriel  fault  is  effectively  de¬ 
coupled  from  the  rest  of  the  fault  system.  If  this  is  true, 
then  from  the  results  indicated  by  the  resolvability  cal¬ 
culations  in  Figure  4.14a,  the  slip  on  this  fault  in  riot 
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detectable  by  the  data.  This  conclusion  Is  not  surprising 
considering  that  we  have  presented  two  separate  dislocation 
models  —  one  without  slip  on  the  San  Gabriel  fault  (Figure 
4.5)  and  one  with  slip  on  the  San  Gabriel  fault  (Figure  4.7) 
—  each  of  which  Is  geophysically  plausible  and  fits  the 
data  to  a  satisfactory  decree. 

4 . 8  Conclusions . 

Several  possible  fault  dislocation  models  have  been 
presented  to  explain  the  observed  vertical  displacement 
which  accompanied  the  1971  San  Fernando,  California  earth¬ 
quake.  The  fault  geometry  assumed  in  calculating  these 
models  was  made  to  conform  to  most  of  the  geological  and 
geophysical  observations  that  were  associated  with  the 
faulting.  In  each  of  the  models  computed  from  the  inver¬ 
sion  scheme,  a  significant  variation  in  slip  along  the 
fault  was  calculated.  The  maximum  displacement  was  com¬ 
puted  to  have  occurred  on  the  shallow  end  of  the  fault 
surface.  The  resolution  of  slip  in  this  shallow  part  was 
very  good,  and  the  spatial  averaging  of  slip  was  relatively 
small.  A  model  was  calculated  in  which  there  was  extensive 
splaying  away  from  the  main  fault  as  the  fault  neared  the 
surface.  The  "ip  model  calculated  for  this  geometry  is 
preferable  to  that  of  a  single  fault  for  the  following  two 
reasons:  1)  The  maximum  displacement  on  the  fault  necessary 

to  fit  the  data  was  reduced,  and  large  variations  in  slip  In 
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the  upper  1  km  of  section  was  reduced;  and  2)  evidence  for 
the  existence  of  and  movement  on  these  splay  faults  was 
observed  in  the  field,  and  their  existence  admits  the 
possibility  of  modeling  the  observed  lateral  offset  in  the 
surface  faulting.  We  caution  here  that  the  exact  number, 
location  and  offset  on  the  splay  elements  should  not  be 
taken  as  absolute.  The  addition  of  these  elements  into 
the  fault  system  was  taken  to  represent  in  some  systematic 
fcShion  the  observable  imbricate  ruptures  which  di^  not 
appear  to  be  so  systematic.  However,  we  think  that  the 
inclusion  of  these  elements  has  led  to  a  better  under¬ 
standing  of  the  faulting  process  near  the  surface. 

Dislocation  models  were  calculated  for  fault  geome¬ 
tries  which  both  included  and  excluded  a  representation  of 
the  San  Gabriel  fault.  It  was  found  that  there  was  very 
little  resolution  of  slip  along  the  San  Gabriel  fault.  If 
slip  did  occur,  this  slip  would  be  normal  fault  movement 
ana  confined  to  the  lower  part  of  the  fault.  These  two 
models  indj  ated  that  between  1^0  and  200  cm  of  slip 
occurred  in  the  area  of  initial  motion.  The  uncertainty  of 
this  slip  in  the  hypocentral  area  as  determined  by  the 
static  data  is  about  100  cm.  All  of  these  calculations 
have  been  based  on  the  initial  assumption  that  the  fault 
sizes  were  fixed.  If  we  have  misjudged  the  length  of  the 
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fault  elements,  the  calculated  slip  would  be  in  slight 
error.  For  distances  away  from  the  fault  elements,  this 
error  is  about  linear  with  length.  That  is,  if  the  lengths 
of  the  fault  elements  in  the  hypocentral  area  are  actually 
10$  smaller  than  the  value  that  was  fixed  in  the  inversion, 
the  calculated  slip  on  the  smaller  fault  element  would  be 
about  10%  larger. 

The  moment  for  this  event  can  now  be  calculated  for 

each  of  the  slip  models.  In  tna  calculations,  the  value  of 

the  shear  modulus  was  taken  from  the  estimates  given  by 

Jungels  (1973,  p.  29).  The  values  of  moment  range  from  1.0 
2  & 

x  10  dyne-cm  for  the  slip  model  presenting  in  Figure  4 . 5 
2  6 

to  2.2  x  10  dyne-cm  for  the  model  given  in  Figure  4.10. 
These  twc  extremes  just  about  span  the  values  obtained  by 
both  static  and  dynamic  observations  given  in  the  litera¬ 
ture  by  various  authors  (Aki,  1971;  Wyss,  1971;  Wyss  and 
Hanks,  1972;  Canitez  and  Toksoz,  1972;  Jungels  and  Frazier, 
1973;  Savage  et  al . ,  1973). 

On  the  basis  of  choosing  the  simplest  possible  model 
necessary  to  fit  the  data,  we  give  as  our  preferred  "best 
fit"  model  the  slip  function  given  in  Figure  4.5  with  the 
addition  of  the  splay  slip  shown  in  Figures  4.10  and  4.11. 
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Appendix  2 

Observed  Gravity  Changes  Associated  with  Deformation 
Acr  -'mpanying  the  1971  San  Fernando  Earthquake 

A2  •  1  Introduction 

In  the  past  there  have  been  a  limited  number  of  stud¬ 
ies  of  gravity  changes  which  accompany  earthquakes.  The 
earliest  of  these  studies  were  in  Japan.  Tsuboi  et  al . 
(1953)  attempted  to  compare  changes  in  gravity  and  precise 
line  leveling  changes  following  the  Nankaido  earthquake  of 
1946.  However,  a  significant  comparison  could  not  be  made 
due  to  the  imprecise  gravity  values.  Extensive  gravity 
surveys  in  the  area  of  the  Matsushiro  earthquake  swarms  of 
1965-1968  were  carried  out  by  Harada  (1968).  A  comparison 
of  this  data  to  the  leveling  resurvey  associated  with  the 
earthquake  swarms  by  Tsubokawa  et  al .  f 19 67)  has  been  re¬ 
ported  by  Kasahara  (1970)  and  Hagiwara  and  Tajirna  (1973). 

In  this  study,  the  relation  between  the  gravity  change  and 
the  vertical  upward  displacement  associated  with  the  tec¬ 
tonic  deformation  of  the  °arthquake  swarms  nearly  coincided 
with  a  Bouguer  gradient  with  a  density  taken  to  be  1.0 
gram/cm  .  The  maximum  displacements  that  were  observed 
where  the  gravity  field  was  being  monitored  were  just  over 
30  cm  and  the  maximum  gravity  changes  were  about  0.0 6  mgal 
with  an  accuracy  of  +  0.02  mgal.  However,  during  the 
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period  of  ground  subsidence  the  changes  in  the  gravity 
field  appeared  to  fall  along  the  free-air  gradient .  As  a 
whole  the  swarm  process  left  no  effective  gravity  change, 
even  though  the  residual  displacement  was  about  20  cm. 

Barnes  (1966)  and  Rice  (1969 )  conducted  remeasure ments 
of  gravity  in  an  area  affected  by  the  deformation  of  the 
Alaska,  196^4,  earthquake.  The  changes  in  gravity  were  com¬ 
pared  with  changes  in  elevation  determined  by  a  post¬ 
earthquake  relevelihg  of  a  1923  first-order  geodetic  line 
coincident  with  the  gravity  traverse.  Because  the  area  in 
which  the  traverse  was  located  had  a  maximum  elevation 
change  of  only  0.6  m,  the  relationship  between  gravity  and 
elevation  changes  could  not  be  determined  accurately.  How¬ 
ever,  the  largest  gravity  differences  indicated  that  the 
slope  of  the  Ag-AE  relationship  was  closer  to  the  Bouguer 
gradient  with  a  density  of  2.67  gram/cm^  than  to  the  free 
air  gradient. 

Hunt  (1970)  has  observed  a  relation  between  changes  in 
the  gravity  field  and  elevation  changes  for  the  1968 
Anangahua  earthquake.  Although  his  data  is  meager,  the 
gradient  of  Ag-AE  relation  between  -0.15  and  -0.20  mgal/m 
is  much  less  than  the  free  air  gradient  -O.309  mgal/m. 

In  each  of  these  cited  studies,  the  conclusions 
reached  by  the  researchers  were  limited  by  either  the  in¬ 
ability  to  accurately  measure  the  vertical  component  of 
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displacement  or  accurately  determine  changes  in  the  gravity 
field  at  points  where  elevation  changes  are  known.  Thus 
in  order  to  establish  a  meaningful  correlation  of  deforma¬ 
tion  and  accompanying  gravity  changes,  the  vertical  dis¬ 
placements  need  to  be  large  and  accurately  known,  and 
reliable  gravity  stations  need  to  have  been  established  in 
the  area  prior  to  the  earthquake.  The  calibration,  stabil¬ 
ity  and  sensitivity  of  gravity  meters  are  sufficient  now 
to  make  gravity  measurements  for  reliable  studies.  Once 
the  correlation  between  AE  and  Ag  has  been  established  then 
the  gravity  meter  could  be  used  as  a  reconnaissance  tool 
to  map  the  extent  of  deformation  rapidly  following  a  large 
earthquake  if  previous  gravity  stations  have  been  estaD- 
llshed . 

Much  of  the  tectonic  deformation  due  to  the  1971, 

San  Fernando,  California,  earthquake  fortunately  occurred 
in  an  area  in  which  there  was  a  dense  leveling  net  and 
enough  reoccupiable  gravity  stations  that  a  study  could  be 
carried  out.  The  gravity  stations  in  the  San  Fernando 
Valley  wer°  installed  as  part  of  a  survey  by  Corbato  (1963) 
using  a  Worden  gravimeter  ( W 88 )  and  very  tight  control  on 
the  observed  gravity  values.  Most  of  the  stations  de¬ 
scribed  by  Corbato  were  actually  installed  in  1958.  These 
stations  were  reoccupied  and  the  results  are  reported  here. 
A  compilation  oi  all  the  observed  gravity  changes  in  the 
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area  affected  by  this  earthquake  is  given  by  Oliver  et  ai. 
(1973)-  This  study  includes  data  in  both  the  San  Fernando 
Valey  and  the  San  Gabriel  Mountains,  both  south  and  east  of 
the  epicenter  and  no.thwest  of  the  earthquake  area. 

A2.2  Gravimeter  Characteristics 

All  gravity  measurements  employed  in  this  study  were 
made  with  a  Worden  Master  (W533)  gravimeter.  The  smail 
dial  constant  of  this  meter  as  determined  from  the  tilt 
tabxe  method  at  the  manufacturer  (Texas  Instruments)  in 
1968  was  0. 3ii6i  mgals  per  large  scale  division  when  the 
instrument  is  ooeroting  at  87°F.  Most  of  the  measurements 
taken  in  this  study  were  taken  when  the  internal  tempera¬ 
ture  of  the  meter  whs  between  70°F  and  72°F.  Since  *be 
respo-- e  of  the  instrument  is  linear.  in  this  temperature 
range,  a  simple  extrapolation  gives  the  dial  constant  at 
72°F  to  be  0.3_>l4o  mgals  per  1  arge  scale  division.  In 
addition,  the  meter  was  checked  on  a  previously  established 
sevet.  station  calibration  loop  between  a  station  locat  a 
at  the  niversity  of  Califox’nia  at  Los  Angeles  (UCI.A)  arid 
Mt.  Wilson.  The  gravity  station  at  the  California  institute 
of  Technology  (Caltech)  is  the  second  station  in  this  Icjp 
denoted  MW2 .  This  calibration  Igod  was  established  by 
arrij'n  and  Cor  ato  C''65)  using  e  LaCoste  and  Romberg 
geodetic  gravimete- .  farther  discussions  and  descriptions 
of  this  loon  can  be  found  in  BieMer  (1964)  and  Robbins 
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e_t  al .  (1973).  The  total  gravity  change  from  MW2  to  the 
station  at  Mt .  Wilson  Is  323*80  mgals.  Successive  checks 
of  the  instrument  over  this  loop  yielded  a  calibration 
value  of  0.3315  +  0.0002  mgals  per  large  scale  division. 
With  careful  measurements  and  averaging  values,  instrument 
readings  at  individual  stations  are  reproducible  to  within 
0.1  large  scale  division,  or  0.03  mgals.  In  the  actual 
survey,  readings  at  individual  stations  were  tied  to  the 
base  station  by  closed  loops,  with  several  readings  being 
made  at  each  of  the  individual  stations. 

The  instruct  drift  due  to  tides,  small  temperature 
changes,  internal  fatigue,  etc.  varied  somewhat  but  never 
seemed  to  exceed  0.15  mgals  per  hour.  The  instrument 
drift  was  linearly  irJ erpoiated  in  time  throughout  the 
closure  of  each  individual  loop. 

A2 . 3  Base  Station  Ties 

All  the  reoccupation  readings  were  based  on  ties  to 
the  base  station  at  Caltech.  The  gravity  station  at 
Caltech,  MW2 ,  used  here  should  not  be  confused  with  another 
close-by  station  described  by  Chapman  (1966,  p.  25,  No. 
309).  A  description  of  the  location  of  this  and  all  other 
stations  referred  to  in  this  chapter  is  given  by  Robbins 
e_t  al .  (1973)’  Individual  reoccupation  stations  were  tied 
to  the  base  station  on  at  least  four  to  five  separate 
occasions,  and  3c  is  thought  that  the  maximum  error  of  each 


tied  station  is  about  +_  0.03  mgals  based  on  the  standard 
deviations  of  the  scatter  in  the  data  about  the  mean  values. 
Two  of  the  reoccupation  stations  tied  to  the  base  station 
were  checked  with  a  LaCos te-Romberg  gravimeter,  with  a 
calibration  constant  of  about  0.1  mgals  per  large  scale 
division,  and  the  differences  in  the  gravity  readings  be¬ 
tween  the  two  stations  and  the  base  station  determined  with 
this  instrument  agreed  to  within  +  0.01  mgal  (Oliver,  per¬ 
son  >mmunication ,  IS 71). 

..  :ce  the  original  gravity  stations  installed  by 
Corbato  were  based  on  direct  ties  to  the  gravity  station  at 
UCLA  and  the  reoccupation  of  these  stations  were  based  on 
ties  to  Caltech,  the  question  arises  as  to  whether  or  not 
there  was  a  relative  change  between  Caltech  and  UCLA  caused 
by  the  earthquake,  or  whether  or  not  the  absolute  gravity 
changed  at  both  of  these  stations.  Oliver  et  al.  (1973) 
reported  that  there  was  indeed  no  gravity  change  at  Caltech 
within  the  reliability  of  5  pre-earthquake  and  6  post¬ 
earthquake  measurements  made  relative  to  UCLA  or  to  a  sta¬ 
tion  at  Menlo  Park,  California.  Supplementary  evidence  for 
this  conclusion  comes  from  the  fact  that  a  continuously 
recording  tidal  gravimeter  operating  at  UCLA  at  this  time 
or  the  earthquake  indicated  no  perceptible  change  in  „he 
gravity  field  at  that  point  relative  to  the  accuracy  of 
that  instrument  (+0.001  mgal).  In  addition,  the  Kt.  Wilson 


calibration  loop  was  rechecked  relative  to  a  distant  cali¬ 
bration  loop,  and  no  changes  in  the  gravity  values  could  be 
observed.  Thus,  we  feel  relatively  comfortable  in  basing 
our  reoccupation  readings  on  ties  to  the  gravity  station 
at  Caltech. 

A2 . 4  Gravity  Station  Reoccupation. 

Figure  A2.1  shows  the  location  and  distribution  of  the 
reoccupled  gravity  stations  in  the  San  Fernando  Valley. 

The  stations  are  identified  with  a  number  corresponding  to 
that  given  by  Corbato  (1963).  The  stations  numbered  1-14 
were  the  sub-base  stations  of  Corbato's  net  and  are  well 
scattered  throughout  the  San  Fernando  Valley.  These  sub¬ 
base  stations  were  the  first  to  be  reoccupled,  since  imme¬ 
diately  following  the  earthquake,  before  any  precision 
leveling  could  be  done,  the  areal  extent  of  the  tectonic 
deformation  was  sought.  It  was  assumed  that  the  deforma¬ 
tion,  if  any,  would  be  reflected  Into  changes  in  gravity. 
The  absolute  gravity  at  these  sub-base  stations  was  more 
accurately  determined  than  at  the  other  stations,  with  the 
values  generally  known  to  within  +  0.03  mgals .  wx  the 
fourteer  stations  originally  comprising  the  sub-base  net, 
eleven  could  be  recovered  in  1971,  alti.cugh  here  was  some 
slight  uncertainty  as  to  the  exact  location  of  three  of  the 
stations.  These  uncertain  stations  are  marked  with  an 
asterisk  in  Table  A2.1.  At  these  uncertain  stations,  an 
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Table  A2.1.  Observed  gravity  and  elevation  changes. 
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attempt  was  made  to  take  the  readings  at  about  the  same 
altitude  as  reported  by  Corbato  ( 1 9 6 3 )  •  The  gravity  values 
for  all  the  stations  for  both  the  pre-earthquake  and  post¬ 
earthquake  measurements  are  listed  here.  The  absolute 
gravity  values  are  the  values  listed  plus  979.0  gals.  The 
values  given  by  Corbato  ( 1 96 3  a  Appendix  B)  are  the  pre¬ 
earthquake  values  rounded  to  the  nearest  0.1  mgal .  However, 
the  values  were  actually  calculated  to  the  nearest  0.01  mgal 
(Corbato,  personal  communication,  1971).  The  uncertainties, 

a,0  ,  are  based  on  standard  deviations  about  the  mean 

gl  g2 

values.  The  uncertainty  of  the  resulting  change  in  gravity, 

0  ,  is  given  by  the  square  root  of  the  sum  of  the  squares 
6 

of  the  standard  deviations  of  the  two  measurements.  The 
change  in  gravity,  Ag,  is  defined  to  be  the  post-earthquake 
measurement  minus  the  pre-earthquake  measurement. 

In  the  last  two  columns,  the  changes  In  elevation, 

AE ,  and  their  respective  standard  deviations  for  the  sta¬ 
tions,  oE,  are  given  when  available.  These  values  were 
obtained  principally  from  the  City  of  Los  Angeles,  Bureau 
of  Engineering,  Van  Nuys  Office.  The  measurements  of  eleva¬ 
tion  both  before  and  after  the  earthquake  are  made  relative 
to  the  mean  tidal  level  in  the  City  of  San  Pedro,  California. 
Since  most  of  the  Corbato  gravity  stations  were  located  on, 
or  very  near,  city  engineering  monuments,  this  data  was 
very  useful.  Supplementary  leveling  information  from 
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Burford  et_  al .  (1971)  was  used.  The  estimated  variance  In 
the  leveling  data  is  from  these  sources.  With  the  exception 
of  two  points  ,  the  changes  in  the  leveling  data  are  rela¬ 
tive  to  the  elevation  determined  in  i960.  These  two  excep¬ 
tions  are  Station  2  and  Station  12  which  are  relative  to 
measurements  in  1970  and  1963  respectively.  Station  12  is 
especially  anomalous  in  that  in  the  period  196 3—19 70  the 
elevation  increased  by  0.06  m,  but  between  1970  and  post¬ 
earthquake  1971 ,  the  elevation  decreased  by  0.06  m. 

From  the  upper  part  of  this  table,  it  is  seen  that  the 

only  significant  changes  in  gravity  occurred  to  the  north 

of  the  observed  surface  ruptures,  with  the  exception  of 

Station  12,  which  is  discussed  below.  Twelve  more  stations 

all  concentrated  above,  or  slightly  below,  the  rupture  zone 

were  then  reoccupied.  The  gravity  values  measured  at  these 

stations  are  listed  in  the  lower  half  of  Table  A2.1. 

A2 . 5  Relation  of  Gravity  Differences  to  Observed 
Deformation . 

The  gravity  stations  listed  in  the  lower  half  of 
Table  A2.1  form  a  profile  from  the  southeast  to  the  north¬ 
west  cutting  across  the  main  surface  ruptures.  The  location 
of  this  profile  in  relation  to  the  surface  deformation  can 
be  seen  in  Figure  A2.H.  Figure  A2.2  shows  the  results  of 
the  reoccupation  of  the  stations  along  this  profile.  The 
top  part  of  the  figure  gives  the  observed  elevation  change 
from  I960  to  1971  showing  the  sharp  scarp  where  the  profile 

8S< 

— . . . .  - - -  -  . . . 


Figure  A2.2.  Gravity  profile  A-A'  showing  the 
relation  of  the  changes  In  surface  elevation  to 
gravity  changes. 
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crosses  the  Sylmar  segment  of  the  surface  breaks.  The 
profile  crosses  just  at  the  western  end  of  the  Tujunga  seg¬ 
ment  of  the  surface  breaks.  The  lower  part  of  the  figure 
indicates  the  corresponding  observed  gravity  change  for 
this  profile.  It  is  clear  that  a  good  relationship  exists 
between  the  observed  elevation  change  and  the  observed 
gravity  change.  Now  if  we  plot  the  significant  changes  in 
gravity  versus  their  respective  changes  in  elevation,  we 
see  that  the  relation  is  quite  linear.  Figure  A2 . 3  shows 
this  plot.  The  line  fit  to  the  data  is  a  least  squares 
relation  given  by 

Ag  =  -  0.215  AE  +  0.026  (A2.1) 

where  AE  is  the  elevation  change  in  meters  and  Ag  is  the 
gravity  change  in  mgals .  The  gradient  of  this  relation 
implies  a  Bouguer  reduction  density  of  2.2  gram/cm  ,  or 
about  the  density  of  the  unrat urated  surficial  layers  in  the 
area,  which  are  late  Cenezoic  sediments.  If  a  Bouguer 
reduction  density  of  0.2  gram/cm^  higher  or  lower  than  this 
value  is  chosen,  the  fit  is  only  slightly  degraded,  with 
the  standard  error  being  increased.  As  an  example,  a 
gradient  of  -  0.224  mgals  per  meter  is  found  when  a  Bouguer 
reduction  density  of  2.0  gram/cm3  is  used.  The  gradient 
for  the  relation  for  the  interval  of  acceptable  Bouguer 


Elevation  change,  m 


Figure  A2.3.  Least-squares  fit  of  elevation  changes 
to  the  gravity  changes.  The  uncertainties  in  .he 
measurements  are  indicated  by  the  error  bars  for 
each  numbered  station.  Points  fall  along  a  gradient 
significantly  less  than  t tie  five-air  gradient. 
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reduction  density  values  is  still  significantly  less  than 
the  free  air  gradient  of  -  0.309  mgals  per  meter. 

If  the  deformation  as  measured  along  this  profile  had 
occurred  as  the  result  of  merely  a  vertical  expansion  over 
some  depth,  then  the  gravity  changes  would  have  fallen  along 
the  free  air  gradient.  However,  the  fact  that  the  values 
fall  along  a  reasonable  Bouguer  reduction  gradient  in  ’.icate 
that  by  some  means  mass  was  added  beneath  the  gravity  sta¬ 
tions.  Within  the  resolving  power  of  the  gravity  measure¬ 
ments,  one  can  conclude  that  the  overriding  thrust  block 
appeared  to  move  as  a  unit.  For  the  regions  of  uplift,  it 
is  clear  that  these  areas  of  negative  gravity  change  are 
now  slightly  heavier  than  before  the  earthquake  due  to  the 
added  mass  from  the  stacking  of  the  surficial  layers  on 
the  pre-earthquake  surface. 

Equation  A2.1  can  be  inverted  to  give  the  estimated 
elevation  change  for  a  known  gravity  change.  This  relation 
is 

AE  =  -il.65  (Ag  +  o  )  .  (A2.2) 

O 

Equation  A2.2  implies  that  if  the  gravity  changes  at  the 
stations  reoccupled  can  be  determined  to  +  0.05  mgals  ,  then 
the  accompanying  change  In  elevation  is  determined  to 


+  0.23  m. 
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The  gravity  change  data  accumulated  in  this  study  has 
been  added  to  that  obtained  by  the  '_>ther  authors  in  Oliver 
et_  al .  (1973)  to  give  a  contour  map  of  the  gravity  changes 
associated  with  this  earthquake  (Figures  A2.4).  Releveling 
data  from  both  that  obtained  from  the  City  of  Los  Angeles 
Engineering  Bureau  and  that  obtained  by  the  U.  S.  Geological 
Survey,  as  reported  in  Burford  et  al^  (1971)  has  been  con¬ 
verted  into  changes  in  gravity  and  included  into  the  gravity 
contours  shown  in  this  figure.  The  gravity  changes  that 
are  obtained  by  conversion  of  leveling  data  are  shown  in 
parentheses  along  side  the  bench  mark  location  and  the  year 
that  the  original  elevation  was  obtained.  Gravity  data 
obtained  directly  are  similarly  shown  without  the  paren¬ 
theses.  Examination  of  the  two  sets  of  data  show  that  they 
are  remarkably  compatible  within  the  uncertainties  of  each 
of  the  data  sets.  Of  particular  interest  is  a  comparison 
1  "  the  gravity  values  from  both  data  sets  in  the  area  Just 
east  of  the  Pacoima  Reservoir.  The  gravity  change  from  the 
converted  elevation  change  (-0.15  mgals)  is  based  on  a 
1929  elevation  determination,  where  the  direct  gravity 
changes  are  based  on  measurements  obtained  in  1964  and 
1970,  yet  the  change  based  on  the  1929  observation  is  about 
the  same  size  as  those  based  on  much  later  observations. 
Savage  et_  al .  (1973)  have  postulated  that  this  particular 
releveling  datum  might  contain  some  pre-earthquake  secular 
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Figure  A2.4.  Contour  map  of  the  gravity  changes 
associated  with  this  earthquake.  This  data 
includes  elevation  changes  that  have  been  converted 
to  gravity  changes  by  equation  (A2.1), 
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uplift,  however,  the  gravity  measurements  here  suggest  that 
most  of  the  uplift  in  this  area  probably  took  place  at  the 
time  of  the  earthquake  and  was  not  distributed  in  the 
years  preceding  the  earthquake. 

Also  indicated  in  Figure  k2.k  are  areas  of  positive 
gravity  changes,  indicating  subsidence.  These  areas  were 
not  delineated  by  the  leveling  surveys  because  of  lack  of 
data.  Although  the  data  indicating  these  positive  changes 
are  less  accurate  than  in  this  study,  the  majority  of  the 
points  defining  the  positive  change  do  have  positive  values 
when  the  maximum  error  is  subtracted.  If  these  gravity 
changes  obey  a  similar  Bouguer  relationship  with  changes  in 
elevation,  then  this  is  an  Indication  that  mass  was  removed 
from  this  area.  No  attempt  has  been  made  to  balance  the 
mass  Involved  because  of  the  poor  nature  of  some  of  the 
data.  However,  the  recognizance  of  this  area  of  subsidence 
shows  the  usefulness  and  facility  of  making  gravity  u„-ure- 
ments  In  tectonic  areas  whether  or  not  they  can  be  done 
in  conjunction  with  leveling  surveys.  It  is  interesting  to 
note  that  this  area  roughly  forms  an  inverted  "U"  and 
corresponds  to  about  the  aftershock  zone  near  the  epicentral 
area  (Allen  et  al^,  1971). 
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A2.6  Bouguer  Residuals. 

Figure  A2.5  shows  the  residual  Bouguer  anomaly  for  the 
profile  A-A'.  This  residual  was  calculated  from  the  for¬ 
mula  , 


=  Ag  -  0.215  AF  .  (A2.3) 

gBA 

The  residuals  are  plotted  in  the  upper  part  of  Figure  A? . 5 . 
Although  the  residuals  are  fairly  small  when  compared  to 
che  estimated  errors  for  each  observation,  the  calculated 
anomalies  are  systematic  in  the  sense  that  the  residuals 
are  negative  south  of  the  surface  ruptures  and  positive  to 
the  north  of  the  surface  ruptures.  Although  barely  signi¬ 
ficant,  these  anomalies  represent  some  rearrangement  of  the 
subsurface  mass  from  1958  to  1971.  One  obvious  rearrange¬ 
ment  of  the  subsurface  mass  is  that  which  can  be  caused  by 
the  raising  or  lowering  of  the  areal  water  table  in  the 
time  interval  between  gravity  measurements.  Oliver  e_t  al . 
(1973)  have  considered  this  in  some  detail.  They  report 
that  for  the  limited  water  table  level  data  that  is  avail¬ 
able  for  this  area,  the  water  table  level  changed  in  the 
short  interval  about  the  time  of  the  earthquake  only  by  the 
amount  the  surface  was  deformed.  Thus  the  entire  sedimen¬ 
tary  unit  sampled  by  the  water  wells  moved  as  a  unit  during 
the  earthquake.  However,  in  the  13-y^ar  interval,  1958  to 


WATER  TABLE 

LEVEL  CHANGE  (1968-1971) 


Figure  A2.5.  Upper — Observed  Bouguer  residuals  for 
the  gravity  profile  A-A'.  Lower— Calculated 
gravity  effects  from  the  fluctuation  in  the  water- 
table  level,  1958-1971,  for  this  profil  . 


1971  >  the  total  water  table  level,  after  seasonal  effects 
are  removed,  was  lowered  In  the  area  north  of  the  surface 
ruptures  and  actually  raised  In  the  area  to  the  south  of 
the  surface  ruptures.  Oliver  (1972,  personal  communication) 
has  calculated  the  gravitational  effects  of  this  net  change 
in  the  water  table  level  and  his  results  for  the  profile 
A-A '  are  plotted  in  the  lower  half  of  Figure  A2.5. 

It  is  seen  here  that  the  water  table  effects  give  rise 
to  about  the  size  gravity  anomaly  as  that  observed,  however, 
it  is  in  the  opposite  direction.  That  is,  the  water  table 
effect  must  be  subtracted  from  the  observed  Bouguer  anomaly 
to  get  the  corrected  total  Bouguer  anomaly.  The  errors  in 
the  gravity  changes  associated  with  the  water  table  fluc¬ 
tuations  are  at  least  as  large  as  the  errors  in  the  upper 
part  of  Figure  A2.5  based  on  the  uncertainty  and  paucity  of 
the  input  data.  It  is  not  wished  to  emphasize  this  anomaly 
too  much,  other  than  perhaps  the  general  form  of  the 
anomaly.  No  attempt  was  made  in  this  study  to  model  the 
anomaly,  since  not  much  is  known  about  the  original  density 
distribution  of  the  subsurface. 

A2.7  Discussion . 

The  reoccupation  of  gravity  stations  after  a  major 
earthquake  is  shown  to  be  a  rapid  and  relatively  inexpen¬ 
sive  method  of  delineating  the  areas  of  deformation  caused 
by  the  earthquake,  provided  a  detailed  gravity  net  is 


established  prior  to  the  occurrence  of  such  an  event. 

The  occurrence  of  the  San  Fernando  earthquake  in  an 
area  where  there  was  such  good  control  on  the  elevation  and 
gravity  changes  has  presented  the  best  opportunity  to 
establish  an  accurate  relation  between  these  two  variables 
for  this  particular  type  of  faulting.  A  gravity  station 
reoccupation  traverse  across  the  zone  of  the  maximum  thrust 
fault  surface  ruptures  yielded  a  propertionality  between 
the  observed  gravity  changes  and  changes  in  elevation.  The 
slope  of  the  relationship  shows  rather  convincingly  that 
the  elevation  changes  and  the  gravity  changes  are  not  re¬ 
lated  simply  by  the  free  air  correction,  but  rather  along 
a  Bouguer  gradient.  The  Bouguer  reduction  density  found  by 
a  least-squares  fit  to  the  data  indicates  that  the  greatest 
portion  of  the  gravity  change  could  be  explaineu  by  the 
stacking  of  the  surficial  layers  on  the  pre-earthquak1': 
surface.  In  retrospect,  this  is  e:actly  what  v.'ould  be  ex¬ 
pected  from  reverse  thrust  faulting.  The  Bougucrs  anomalie 
derived  for  this  profile,  although  statistically  barely 
significant,  are  systematic  In  that  they  are  generally 
positive  to  the  north  of  the  rupture  zone  and  negative  to 
the  south  of  the  rupture  zone. 


